{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Refraction seismics - fitting a layered earth model\n",
    "O. Kaufmann - Nov. 2017 - GPL v.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from collections import OrderedDict\n",
    "from scipy.optimize import least_squares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "picked_times = ''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_span = [-1., 50.]\n",
    "sources = [0., 23.5, 49.]\n",
    "receivers = np.arange(1., 49.)\n",
    "v = [254., 1487]\n",
    "z = [0, -3.7]\n",
    "theta = [0.0, -0.05]\n",
    "td = [0.003, -0.002, 0.001]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class layered_earth_model(object):\n",
    "    def __init__(self, v, z, theta, td):\n",
    "        '''\n",
    "        v : a list of compressive or shear waves velocities in the layers\n",
    "        z : a list of the elevation of the top of the layers at x=0 \n",
    "        theta : a list of dipping angles (in radians counterclockwise) of the top of the layers\n",
    "        '''\n",
    "        self.v = v\n",
    "        self.z = z\n",
    "        self.theta = theta\n",
    "        self.td = td  # trigger delays\n",
    "    \n",
    "    @property\n",
    "    def v(self):\n",
    "        \"\"\"Gets and sets the velocities in the model\n",
    "\n",
    "        \"\"\"\n",
    "        return self.__v\n",
    "    \n",
    "    @v.setter\n",
    "    def v(self, val):\n",
    "        try:\n",
    "            self.__v = list(val)\n",
    "        except:\n",
    "            print('could not set velocities')\n",
    "            \n",
    "    @property\n",
    "    def z(self):\n",
    "        \"\"\"Gets and sets the elevations in the model\n",
    "\n",
    "        \"\"\"\n",
    "        return self.__z\n",
    "    \n",
    "    @z.setter\n",
    "    def z(self, val):\n",
    "        try:\n",
    "            self.__z = list(val)\n",
    "        except:\n",
    "            print('could not set elevations')\n",
    "    \n",
    "    @property\n",
    "    def theta(self):\n",
    "        \"\"\"Gets and sets the dipping angles in the model\n",
    "\n",
    "        \"\"\"\n",
    "        return self.__theta\n",
    "    \n",
    "    @theta.setter\n",
    "    def theta(self, val):\n",
    "        try:\n",
    "            self.__theta = val\n",
    "        except:\n",
    "            print('could not set dipping angles')\n",
    "    \n",
    "    def layers_count(self):\n",
    "        ''' Get the number of layers in the model\n",
    "        \n",
    "        :return: Number of layers in the model\n",
    "        :rtype: float\n",
    "        '''\n",
    "        return len(self.v)\n",
    "    \n",
    "    def plot_model(self, x_span=(0., 1.), filled= True, *args, **kwargs):\n",
    "        x_span.sort()\n",
    "        z = [[self.z[i] + x_span[0] * np.tan(self.theta[i]), self.z[i] + x_span[1] * np.tan(self.theta[i])] \n",
    "             for i in range(self.layers_count())]\n",
    "        for i in range(self.layers_count()):\n",
    "            plt.plot(x_span, z[i], '-k', *args, **kwargs)\n",
    "        if filled:\n",
    "            for i in range(self.layers_count()-1):\n",
    "                label = 'V%1d : %.0f m/s' %(i, self.v[i])\n",
    "                plt.fill_between(x_span, z[i-1], z[i], label=label)\n",
    "            label = 'V%1d : %.0f m/s' %(self.layers_count(), self.v[-1])\n",
    "            plt.fill_between(x_span, z[-1], [1.2*min(z[-1]), 1.2*min(z[-1])], label=label)\n",
    "            \n",
    "    def plot_layout(self, sources, receivers):\n",
    "        plt.plot(sources, np.zeros(len(sources)), '*r')\n",
    "        plt.plot(receivers, np.zeros(len(receivers)), 'vk')\n",
    "        \n",
    "    def plot_response(self, sources, receivers, times, *args, **kwargs):\n",
    "        plt.plot(sources, np.zeros(len(sources)), '*r')\n",
    "        for i in range(len(sources)):\n",
    "            plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
    "            \n",
    "    def two_layers_forward_refraction(self, sources, receivers):\n",
    "        ic = np.arcsin(self.v[0]/self.v[1])\n",
    "        x_span = [min(min(sources), min(receivers)), max(max(sources), max(receivers))]\n",
    "        times = np.array([np.zeros(len(receivers)) for i in range(len(sources))])\n",
    "        for j in range(len(sources)):\n",
    "            h = (self.z[0]-self.z[1] - sources[j] * np.tan(self.theta[1]-self.theta[0])) * np.cos(self.theta[1]-self.theta[0]) # Thickness of layer one at the source point\n",
    "            t_direct = [abs(receivers[i]-sources[j])/self.v[0] for i in range(len(receivers))]\n",
    "            t_refracted = [(abs(receivers[i]-sources[j])/self.v[0])*np.sin(ic-np.sign(receivers[i]-sources[j])*(self.theta[1]-self.theta[0])) + (2 * h * np.cos(ic))/self.v[0] for i in range(len(receivers))]\n",
    "            stacked = np.dstack((t_direct, t_refracted))\n",
    "            times[j] = stacked.min(2).flatten() + self.td[j]\n",
    "        return times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7faa341586a0>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsEAAAD4CAYAAAANWzs4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAoIElEQVR4nO3deXyU5b338e9vJhskE8gqCCLuu4JGUayPYFGg1KqttdJFPVqpnrZHW63HrdZq7WYt+Dyn1WLrsa0eOdZWsW6guFWtC1iqWMQFwZmwhSWZsISQ5Hr+mEkyk8xMAjNJJnN/3q8XLzL3XHPf13jD+PXnb67LnHMCAAAAvMQ30BMAAAAA+hshGAAAAJ5DCAYAAIDnEIIBAADgOYRgAAAAeE7eQFy0srLSjR07diAuDQAAAA9ZsmTJRudcVdfjAxKCx44dq8WLFw/EpQEAAOAhZrY60XHaIQAAAOA5hGAAAAB4DiEYAAAAnkMIBgAAgOcQggEAAOA5GQnBZjbNzFaY2Ydmdm0mzjng1q6VTj1VWrcu5bDx48fLzLr9Gj9+POMZn3J8Ns2F8QMwvstnzIDPpx/HZ9NcGM/4dMZn01z6Y3yHXmakbJd2CDYzv6RfSZou6XBJM83s8HTPO+BuvVV6+WXplltSDjvppJNUUFAQd6ygoEATJ05kPONTjs+muTB+AMZ3+YwZ8Pn04/hsmgvjGZ/O+GyaS3+M79DLjJTtzDmX3gnMTpJ0s3NuavTxdZLknPtJstfU1NS4/l4n+Morr9TSpUt7HLfgb39TYVtbt+M7fT5NPeWU7sd37tQbb7yhtpjX+Hw+TZgwodsfLMYzPlvnwvj+G7/SOe2b4HN3lZkOMMv6+fNnn/GMz8659Mf43c1IscaNG6c5c+akHNNXzGyJc66m6/FMbJYxSlIw5nFI0oQEE5glaZYkjRkzJgOX3X1r6neofvuulGOO2/dw3VxXq89sq9dQ57TdTE8UD9cPq0Zpw5pwwtcUlJSpKbxZkpNkKgiU6cONTZKaGM/4lOOzaS6M75/xM8pHJv2MKdi0Nuvnz599xjM+e+fS1+P3JCO1O/SI1pTPD4R+2zHOOTdX0lwpUgnur+u2mzNnjgoe/qceWhxKOW6rpKYFv1LR0qfV5M9XUWuLmg48UVunflNDk7ymYOtmrfnN1+VammV5+ar86hz5S8qSXoPxjM/GuTC+f8ZvLSlL+hlTMQjmz599xjM+e+fS1+P3JCO1u+na03oY0f8y8cW4Wkn7xDweHT02aFVur9f946frnAvu0P3jp6tqe33K8Xkl5So+8tOSTMVHTUn5h43xjM/WuTC+/8Yn+4wZLPPPxPhsmgvjGZ/O+GyaS3+M392MlM0yUQl+U9JBZrafIuH3fElfzsB5B8xl59zQ8fNNZ/x7r14z7OSZ2rUpqOETZzKe8bs1Ppvmwvj+GZ/qM2YwzD9T47NpLoxnfDrjs2kufT1+TzJStkr7i3GSZGafkTRHkl/Svc6521KNH4gvxknSNb1ohwAAAEBmvXrtadp7+JABuXZffjFOzrknJT2ZiXMBAAAAfY0d4wAAAOA5hGAAAAB4DiEYAAAAnkMIBgAAgOcQggEAAOA5hGAAAAB4DiEYAAAAnkMIBgAAgOcQggEAAOA5hGAAAAB4DiEYAAAAnkMIBgAAgOcQggEAAOA5hGAAAAB4DiEYAAAAnkMIBgAAgOcQggEAAOA5hGAAAAB4Tloh2My+aGbvmlmbmdVkalIAAABAX0q3ErxM0uclvZSBuQAAAAD9Ii+dFzvnlkuSmWVmNgAAAEA/6LeeYDObZWaLzWxxXV1df10WAAAA6KbHSrCZPStpRIKnbnDOze/thZxzcyXNlaSamhrX6xkCAAAAGdZjCHbOTemPiQAAAAD9hSXSAAAA4DnpLpF2jpmFJJ0k6QkzW5CZaQEAAAB9J93VIR6R9EiG5gIAAAD0C9ohAAAA4DmEYAAAAHgOIRgAAACeQwgGAACA5xCCAQAA4DmEYAAAAHgOIRgAAACeQwgGAACA5xCCAQAA4DmEYAAAAHgOIRgAAACeQwgGAACA5xCCAQAA4DmEYAAAAHgOIRgAAACeQwgGAACA5xCCAQAA4DlphWAzu93M3jOzt83sETMbnqF5AQAAAH0m3UrwM5KOdM4dLel9SdelPyUAAACgb6UVgp1zC51zLdGHr0kanf6UAAAAgL6VyZ7giyU9lexJM5tlZovNbHFdXV0GLwsAAADsnryeBpjZs5JGJHjqBufc/OiYGyS1SHog2Xmcc3MlzZWkmpoat0ezBQAAADKgxxDsnJuS6nkzu0jSZyV92jlHuAUAAEDW6zEEp2Jm0yRdI+lU59z2zEwJAAAA6Fvp9gT/l6SApGfMbKmZ3Z2BOQEAAAB9Kq1KsHPuwExNBAAAAOgv7BgHAAAAzyEEAwAAwHMIwQAAAPAcQjAAAAA8hxAMAAAAzyEEAwAAwHMIwQAAAPAcQjAAAAA8hxAMAAAAzyEEAwAAwHMIwQAAAPAcQjAAAAA8hxAMAAAAzyEEAwAAwHMIwQAAAPAcQjAAAAA8hxAMAAAAzyEEAwAAwHPSCsFmdquZvW1mS81soZntnamJAQAAAH0l3Urw7c65o51z4yQ9Lumm9KcEAAAA9K20QrBzLhzzsFiSS286AAAAQN/LS/cEZnabpAskNUianGLcLEmzJGnMmDHpXhYAAADYYz1Wgs3sWTNbluDXWZLknLvBObePpAckfSvZeZxzc51zNc65mqqqqsy9AwAAAGA39VgJds5N6eW5HpD0pKQfpDUjAAAAoI+luzrEQTEPz5L0XnrTAQAAAPpeuj3BPzWzQyS1SVot6bL0pwQAAAD0rbRCsHPuC5maCAAAANBf2DEOAAAAnkMIBgAAgOcQggEAAOA5hGAAAAB4DiEYAAAAnkMIBgAAgOcQggEAAOA5hGAAAAB4DiEYAAAAnkMIBgAAgOektW3yYLJkyRI9d99shbcXyR+oVF5ppfyBSvlLymU+/0BPDwAAAP3IMyH4/fff13uL/qRtTc1xx81MRaVlUkmVrKRK/kCF8kqrIkE50B6UywjKAAAAOcQzIXjmzJk6f8jf1PDa/Qo2tCkUblMw7KK/NyoUblBw80oFP25VY3Nb3GvN51NRYHhnUC7tDMh5gUr5SyvlLyYoAwAADBaeCcFSpOo7vMg0vMivo/ZKHFidc6pvUjQctykUdpHQ3NioYEODQimDcplUUikLVEcqyoHYynIFQRkAACBLeCoE94aZqWyIVDak56AcCcltCja0V5TDCoXrFdy4UsGPWtW4K0lQDkQqyh3tFjGVZX/xcIIyAABAHyME74HYoHx0iqC8pb2i3BCtKIfbFAqHFQzXK7TxIwU/alW4S1D2+XwqLC2PVJTbg3JpZXyPMkEZAAAgLYTgPmJmKh8ilfciKHfvUW5QKLxFwWhQbkoalKtkJZUd7Rbt7Rf+0qpIUDZWwAMAAEgkIyHYzK6S9AtJVc65jZk4pxfEBuVjRiQPypt3uM6A3NAlKNd9qFCSoFxUWi5XUiULVHX/Il97RZmgDAAAPCjtEGxm+0g6Q9In6U8HXZmZKoaaKoZK41IE5U07XOeX+DpaLxoUDG9RaMOHCn2YqqJcLYtpt8jraL+okq94GEEZAADknExUgmdLukbS/AycC3vAzFQ51FTZi6DcWUnuUlHe8IFCH7RoZ4uLe53P71dRabnaiqvkC1TFBOTOwOwbSlAGAACDS1oh2MzOklTrnPunmfU0dpakWZI0ZsyYdC6LPRAblMePTB6UN253MZXk9uXh6hVs2KxQkqDs9+epsLRcbSWVkaAc13pRpbxABUEZAABklR5DsJk9K2lEgqdukHS9Iq0QPXLOzZU0V5JqampcD8MxAMxMVcWmquKeg3L3HuUtCoU3KbjuA4Xeb1FzkqAc6VGu7LKGciQw+4aWEpQBAEC/6DEEO+emJDpuZkdJ2k9SexV4tKS3zOwE59y6jM4SWSM2KB+bIijXbU/Uo7xFwfAmhdb5egzKvkBVXH9yZ1Aepp7+rwMAAEBP9rgdwjn3jqTq9sdmtkpSDatDwMxUXWyqThGU29oryt16lDdHWi/Wva/Qihbtak3WelHdZaORzsqyb0gpQRkAAKTEOsEYEL6YoHzc3smDct22RK0XmxUKb1Zwrak2UVDOy1NhaYXaiqsiayZ3rKHcWVkmKAMA4G0ZC8HOubGZOhcgRYLyXiWmvUqkmhRBecO2RK0XmyKtF2tWqDacIii3V5Rj1lBu/2IfQRkAgNxFJRiDms9MI0pMI3oRlDsrye1heVMkLK8x1Ta0qKUtPijn5eWroLRCbSVVnUG5y+58BGUAAAYnQjByXmxQPn5U8qC8fmvs8nDR9ovGjQqGNyrU26Bc2n13Pl9RgKAMAECWIQQDigTlkQHTyEDyoNzaFq0od/Qot1eUN0bCcq2pdnmLWrsG5fyCSFAuroz2KHffnY+gDABA/yIEA73k93UG5RNSBOX1cT3K7Stf1CkYrlMo9J5qwymCckfrRWSTkfbNRvyBKvmKSgjKAABkCCEYyCC/z7R3wLR3L4Jy/Bf5XCQkh+sUDJnWJArKBQUqCHR+ma/rF/n8pVXyFRYTlAEA6AVCMNDPOoOyTxOSjGltc1q3tcv21e1BubFOwaC0JtyiLjk5EpRLKyPLw3WtJkdbMQjKAAAQgoGs5PeZRpWaRpUm30a6JRqU49dQdgqFNygY3qBQkqCcX1Co/PbWi5Kqjr7k2MqyEZQBADmOEAwMUnk+0+hS0+hSn04cnXhMe1Du3nqxIRKWg9LapEE5WlGOabkgKAMAcgUhGMhhsUE5mZY2p7WNXbavbnAKNa5XMLxeoU+SBOXCQuUHKju+zJcXs3V1++9WMJSgDADISoRgwOPyfKZ9hpn2GebTSUnG7GpN1qO8PhKWV0trG1vkugXlokjrRUePcvyX+fJKIz3KAAD0N0IwgB7l+zuDcjK7Wp3WdutRblMoHK0opwrKsatexPYoR3fp8xUO7eN3CADwGkIwgIzI95vGDDONGeaT9kk8pj0ox/cot8VVlNclrShHepTb2y3iV7+oJCgDAHYLIRhAv4kLyknsanVaE9Oj3LnpyDqFGtcp+LG0fmv3oFxQNER50Ypyx9rJXSrLBGUAQDtCMICsku837TvctO/w5EG5OTYod7RetCoUXqdg4zqFPo5UlLuKBOXIl/nyohuM5HVUlSP9yr6CIX359gAAWYIQDGDQKfCbxg43je1FUI5vvWhVMLxWoca1kYpysqBcWqW24spo60X8F/n8AYIyAOQCQjCAnNTboFwb7ro8XItCjWsUDK9RaGWSoDxkaHxFuWMN5c7Ksq+gqC/fHgAgTYRgAJ5V4DftV2baryx5UN7ZEq0oxy0PtysSkhvXKLhB2rC1p6Acu4ZyZ2WZoAwAA4cQDAApFOb1LijXdutR3qVQeI2CjWsUShKUC4cUyx+okIsuD5dodz5fPkEZAPpCWiHYzG6WdKmkuuih651zT6Y7KQAYTArzTPuXmfZPEZSbWhL1KDfHVZTrkgXl0kq56BbWiXbnIygDwO7LRCV4tnPuFxk4DwDkrKJeBuXOHuX29otmhRprFQzXKpQsKA8tlj9QKRe3hXVsZZmgDABd0Q4BAFmiKM90QLnpgPLUQTnU8SW+9qpys0LhWgUbaxVaJ23cliwoV8UE5cqYynK0Rzm/sC/fHgBklUyE4G+Z2QWSFku6yjm3JdEgM5slaZYkjRkzJgOXBQDvKcozHVhuOjBFUN6xK9Kj3LnRiIu2XoQUagwpuM5p07bWbq8rHFoSX1GOtlvkBaKtGCUVBGUAOaPHEGxmz0oakeCpGyTdJelWSS76+x2SLk50HufcXElzJammpsYlGgMASN+Q/N4F5VA4tj/ZKdiwMxKSwyGF1kmbElWUiwNdWi+678xneQV9+fYAICN6DMHOuSm9OZGZ3SPp8bRnBADoc0PyTQdV+HVQRfIx23c51catoewUCjcpGA4qtDWk4Jo2bd6eoKIcF5SjayfH7c5XQVAGMODSXR1ipHNubfThOZKWpT8lAEA2GNrLoByKW0PZKRhuUigcVLAxpFCSoFxUHJAFqqSO1ov4/uRIRTm/D98dAK9Ltyf452Y2TpF2iFWSvpHuhAAAg8fQfNPBFX4dnCIob2t2qm2MXUO5vaL8iUKNQQXXtGlLkqDsi36Zrz0kx+3OV1JBUAawx9IKwc65r2VqIgCA3FRc0Lug3L1HuUmhxk8UbAwqlCwol5TKSiqlQHVHu0Xc7nwEZQBJsEQaAGDAFReYDqn065DK5GO2NidqvdihUPiTSJ9yqE31O5IE5WjrRdcv8hGUAe8iBAMABoWSAtOhlX4d2ougHN96sUPB8KpIWE4alIfJApWdQTlmC+tIC0aFzE9QBnIJIRgAkDN6E5QbdyZqvdiuUOOqSJ9yyqBcJZV0WUO5vf2ipJygDAwihGAAgKcECk2HVfl1WFXyMe1BOX55uO0Khj9WKLxawWCb6pvig7KZqbCkVFZSJQWq4la6iFSW24My/+oFsgF/EwEA6KI3QTm8M1GP8naFwh8rWL9aoU8SB+WikmFSSZUsWkWOX0u5kqAM9BP+lgEAsAdKC02HV/l1eC+CcmePcptC4W0KhhsVql+l4Oo21e9MEpQDVbLoOspxaymXVspfTFAG0sXfIAAA+khvg3JnJbktGpq3KdTYqGD9KgVXt6p+Z1v8i8w0JDBcKqnsCMqdK15UKa+0Qv6SCpnP36fvDxjMCMEAAAyg0kLTEdV+HVGd+HnnnMI7FfMlvmhQDm9VKBxWcMsqBVe1aktzfFA2MxUFhkdaL0oq47au7gjNJeUEZXgWIRgAgCxmZhpWJA0r8uvIFEG5YaeiVeTYqvJWBcNhhbZ83HNQbg/HXddSJigjRxGCAQAY5MxMw4uk4UV+HVmdOLC2B+XOSrKLryhv/ljBj1vVmCwoB6ojFeXYravbK8slZQRlDDqEYAAAPCA2KB+1V/KgXN+kmC/xRdsvGrcq2BBWaPPKxEHZ54tpveiydXX7piMEZWQZQjAAAJAUCcplQ6SyIT0H5c4v8bVXlBsVCjcouGmlgitb1birh6Acsytf+1rK/mKCMvoPIRgAAPRabFA+OkVQ3tKUqEe5UcFwg0Ipg3JZZNWLQHXnrnyx21gXDycoIyMIwQAAIKPMTOVDpPJeBOXuPcphhcL1Cm5cqeBH3YOyz+dTYWl5x/JwHVXkAEEZu4cQDAAA+l1sUD5mRPKgvHmH67J9dVukmhzeomDdRwp91KqmpEE5sjxcpD85WlWO7s7nLx4uM19/vFVkKUIwAADISmamiqGmiqHSuBRBedMO12X76rZIf3J4i0J1HyYNykWl5XLR5eG69Se3V5QJyjmLEAwAAAYtM1PlUFNlL4JyZyW5S0V5w4cKfdCinS0u7nWdFeVqWYI1lPMCVfIVDyMoD1Jph2Az+7akb0pqlfSEc+6atGcFAACQIbFBefzI5EF543YXU0luXx6uQcGGLQpt+CBxUPb7VVRarrbiKvkC0bWTY9ZQzgtUEpSzVFoh2MwmSzpL0jHOuZ1mlmQvGwAAgOxlZqoqNlUV9xyUu/co1ysU3qzg+khQbk4WlEuq5Cupiqkkd1aWfUMJyv0t3Urw5ZJ+6pzbKUnOuQ3pTwkAACD7xAblY1ME5brtiXqU6xUMb1YoSVD2+/NUWFqutpLKaEW5++58BOXMSjcEHyzpFDO7TVKTpKudc2+mPy0AAIDBx8xUXWyqThGU29oryt16lDcr2LBJoXUfKLSiRbtaEwflyJf5KjtWu4jdnS8SlK0/3uqg12MINrNnJY1I8NQN0deXSzpR0vGSHjKz/Z1zrutgM5slaZYkjRkzJp05AwAADFq+mKB83N7Jg3LdtkStF5sjrRdr31dtiqDcVlIlf9xGI52B2TeklKCsXoRg59yUZM+Z2eWS/hINvW+YWZukSkl1Cc4zV9JcSaqpqekWkgEAABDhM9NeJaa9SqSaFEF5w7ZErRebI2E5WVDOy1NhaYXaiqsiayZ3rKHcufKFF4Jyuu0Qj0qaLOl5MztYUoGkjelOCgAAAKn5zDSixDSiF0G5s5LcHpY3KRTepOCaFaptaFFLW5KgXFIdWTM5Zg3l9n7lwR6U0w3B90q618yWSWqWdGGiVggAAAD0v9igfPyo5EF5/dbY5eGi7ReNmyJheY0lDMp5efkqKK2Itl60r3QRvztfNgfltEKwc65Z0lczNBcAAAD0M5+ZRgZMIwPJg3JrW7Si3NGj3F5R3qhQ40YFa6XacKtakwTlS/81QU889mg/vJveY8c4AAAApOT3dQblE1IE5fVxPcrtK19sVFFJfj/PuGeEYAAAAKTN7zPtHTDtnSgof2f2wEwqBVZcBgAAgOcQggEAAOA5hGAAAAB4Ttb0BO/atUuhUEhNTU19d5HRX5aqP9d3589ZTkUNKzX6rZ8pv7l+oCcDAACQtqwJwaFQSIFAQGPHju279eTqV0vbN/fNuXOYc06btpUrpP/Ufq9dN9DTAQAASFvWtEM0NTWpoqIiaxdU9jIzU0VxnpqG7T/QUwEAAMiIrAnBkgjAWSxyb7g/AAAgN2RVCAYAAAD6Q9b0BHc19tonMnq+VT+dkfL5yefO0rXfukhTJ03sODbnnge04qPVuuun1+v3D/1VP7rzt5KkG6/4ui4878xeX/uBvzypn/36PjknBYqH6q6fXK9jjjhYkjR2wgwFSorl9/mUl+fX4qceiHvtHXf/UVffOlt17yxSZXlZr6/Z1WtL3tbv5s3XPbd/f4/PAQAAkCuyNgT3t5lnT9W8+QviQvC8+Qv08xuv0OYtDfrh7Lla/OT9MjMdN/0r+twZp6pseGmvzr3fPqP04sO/VdnwUj313Cua9Z8/0uuP/6Hj+ef/9JuEATdYu04LX/q7xowakfb7e+r5VzQt5r0BAAB4Ge0QUefOmKInFr2s5uZdkqRVwTVas36jTplwrBa8+HedfsoElZcNU9nwUp1+ygQ9/cKrvT73xOOP6QjMJx57lEJr1/fqdd+5+Q79/IYrk/ZK3/e/j+nsi7+r08+/XGMnzNB//fc8/fI392v8GTN14mcv0OYtDR1jF738pqaccoLeXfGRTpjxNY07/XwdPeU8fbDyk16/DwAAgFxBCI4qLxumE8Ydoaeef0VSpAp83pmny8xUu26D9tm7sxo7euReql23ods5brr9Lj228MWU1/ndvEc1ffLJHY/NTGfM/KaOm/Zlzb3/zx3H5y94QaNGVne0TSSzbMWH+stvf6E3n7xfN/zs1xo6pEj/WPigTjruaP3h4cclSRs3b1F+Xp6GlQZ09x8f1hWXzNTSZ+Zp8ZMPaPTI6p7/4QAAAOQY2iFizDx7mubNX6Czpk7SvPkL9Ls7btqt19/yvctTPv/8K2/qdw8+qpcfubfj2MuP3KtRI6u1YeNmnX7+5Tr0wLGqOeZw/fj/3auF//OrHq85eeLxCpQUK1BSrGGBEp15+v+RJB112IF6+18fSJIWvviazjj1REnSSccdrdv+7+8UWrtBn59+mg7af8xuvUcAAIBcQCU4xllTJ2nRy2/orXeWa/uOJh139OGSpFEjqhVcs65jXGjteo0asXsV1Lf/9b6+/r1bNf/e2aooH95xfFS0EltdWa5zpk/WG0vf1UerQvr4k1odc/r5GjthhkJrN+jYqV/Rug0bu523sCC/42efz1RYGHnsM59aWlslSU8994qmTY70A3/5nOl67L9na0hRoT7ztW/ruZff2K33AQAAkAsIwTFKiodq8sQaXfzdH2rm2VM7jk899SQtfOk1bakPa0t9WAtfek1TTz2p1+f9pHatPn/p1frjnbfq4AP27Ti+bfsONW7d1vHzwhdf05GHHKCjDjtIG95epFWvP6FVrz+h0SOr9daCBzSiunK335NzTm8v/0DjjjhEkrRydUj77zta/3HJTJ01dZLeXv7Bbp8TAABgsMvadoieljTrKzPPnqZzLrlK8+76Scex8rJh+v6VX9fxM74qSbrpO5eqvGxYt9fedPtdqjnmcH3ujFPjjt8y+x5t2tKgf78+cs72pdDW123SOZdcJUlqaW3Vl8+epmkx/cKZsOTt5Rp/5KEdX6576K/P6I9/fkL5eXkaUV2h6799cUavBwAAMBiYc27PX2z2v5IOiT4cLqneOTeup9fV1NS4xYsXxx1bvny5DjvssD2eS6/Ur5a2b+7ba2SZH835rQ7cbx+df9bUngf3YPnqDTpswXkZmBUAAPCU77wrDRs9IJc2syXOuZqux9OqBDvnvhRzgTskNaQYjgFw45VfH+gpAAAAZJ2MtENY5P+1nyfptEycDwAAAOhLmfpi3CmS1jvnkn7LysxmmdliM1tcV1eXocsCAAAAu6/HSrCZPSsp0b69Nzjn5kd/ninpwVTncc7NlTRXivQE7+Y8AQAAgIzpMQQ756aket7M8iR9XtJxmZoUAAAA0Jcy0Q4xRdJ7zrlQBs4FAAAA9LlMfDHufPXQCrFHbu6+Dm9650u9cMXkc2fp2m9dpKmTJnYcm3PPA1rx0Wp946tf0OXX/Vjhrdvk9/t0w7cv0Zd2Y8mxTZvrde6sa/TmP9/VReedqf+67dpuYz530ZVa+Umtlj33J0nS0mUrdNm1t6lpZ7Py8vz69Y+v0wnjj9Ttd/1eD/zlKUmRtYWXf/Cx6t5elHDd4t46btqX9er8+1RYWLDH5wAAABhM0g7BzrmLMjCPATfz7KmaN39BXAieN3+Bfn7jFRo6pEh/uPNWHbT/GK1ZV6fjpn9FUydN1PBhgV6du6ioULdec7mWvfeRlq34sNvzf3lykUqKh8Ydu+a2O/WD735D0087WU8uelnX3HanXnj4Hn3v8gv1vcsvlCT9deGLmn3PA2kF4I8/qdWoEdUEYAAA4Clsmxx17owpemLRy2pu3iVJWhVcozXrN+qUCcfq4AP21UH7j5Ek7T2iStUVZarbtKXX5y4eOkSfOmG8ihIEza3btuuXcx/QjVfEr+drJoUbt0qSGhq3au+9qrq99sH5CzTz7GkJr1ly0Mn63q2zdcTkczXlS5fpjX8s06RzL9X+J52pxxa+2DHu6edf1bRJE9Xa2qqLrvyBjjztizrq0+dp9tz7e/3+AAAABhtCcFR52TCdMO4IPfX8K5IiVeDzzjy9Y7vhdm/8Y5mad+3SAWO773py9x8e1t1/eHi3rvv9n/9aV33jqxo6pCju+JwfXq3v/ehO7VMzXVffOls/ue5bcc9v37FDT7/wqr7wmU8nPO+27Tt02skn6N3nH1agpFg3/vzXeubBX+uR396hm26/q2Pc0y+8qmmTJ2rpuytUu26Dlj33J72z6CH925fO2q33AQAAMJgQgmPMPHua5s1fICkSgmeeHd/3u3Z9nb72H9/Xf//yZvl83f/RXXbBubrsgnN7fb2ly1boo9UhnTO9+x4jd/3hYc2++SoFFz+l2T+4SpdcdUvc839d+JJOrjkmaStEQUG+pk2OtHYcdeiBOvXEY5Wfn6+jDjtQq0JrJEnNzbsUWrte++87WvuPGa2Vn9Tq2zf+TE8//4pKA8W9fh8AAACDDSE4xllTJ2nRy2/orXeWa/uOJh139OEdz4Ubt2rGBVfotv/8pk487uiMXO/vS97W4rf/pbETZuhTZ1+s91eu1qRzL5Uk/f5Pj+vzn4mE4y+eebreWPpu3GvnPbYwaSuEJOXn5XVUsX0+6+j59fl8amlplST97fW39KkTxkuSyoaX6p/PzNOkk2p09x//rK9ffUviEwMAAOQAQnCMkuKhmjyxRhd/94dxVeDm5l0655KrdMG5M3TuZ1Mum7xbLr/wi1rz1kKtev0JvfzovTp4/331wsP3SJL23qtSL/59iSTpuZff0EH77dPxuoZwo158bYnOmjopres//cKrmh6tFm/cvEVtbW36woxP60fX/Lveeue9tM4NAACQzTKxRFrf6GFJsz1ifsmX+i3PPGeGzrn4Ss27+/aOsQ89/rReev0f2lQf1n1/elySdN+cWzXuyEPjXnv37x+SJF124Xndzjv2+GkKb92q5uZdenTBC1r44G90+CEHdA7w+aO/R655zy9u1hXf/5laWltVVFigubff3PHcIwte0hmnTlRxSQ+rU7S/V/NFfsW+d1+eXvj7W7rlmm9LvjzVrt+sf7vy+2pzkc38fnL9Fd3/WZlfGlqR+poAAABdWfbVXc25/t/BuKamxi1evDju2PLly3XYYYf1+1y8KhQK6dJLL9VTTz3V69dwjwAAwGBjZkucczVdj2dfLEe/GD169G4FYAAAgFxCCAYAAIDnZFUIHojWDPQO9wYAAOSSrAnBRUVF2rRpE2ErCznntGnTJhUVFfU8GAAAYBDImtUhRo8erVAopLq6uoGeChIoKirS6NHdd8kDAAAYjLImBOfn52u//fYb6GkAAADAA7KmHQIAAADoL4RgAAAAeA4hGAAAAJ4zIDvGmVmdpNX9fuHdVylp40BPAn2G+5u7uLe5i3ub27i/uWsg7+2+zrmqrgcHJAQPFma2ONE2e8gN3N/cxb3NXdzb3Mb9zV3ZeG9phwAAAIDnEIIBAADgOYTg1OYO9ATQp7i/uYt7m7u4t7mN+5u7su7e0hMMAAAAz6ESDAAAAM8hBAMAAMBzCMFJmNk0M1thZh+a2bUDPR/sOTO718w2mNmymGPlZvaMmX0Q/b1sIOeIPWNm+5jZ82b2LzN718yuiB7n/uYAMysyszfM7J/R+/vD6PH9zOz16Ofz/5pZwUDPFXvGzPxm9g8zezz6mHubI8xslZm9Y2ZLzWxx9FhWfTYTghMwM7+kX0maLulwSTPN7PCBnRXScJ+kaV2OXStpkXPuIEmLoo8x+LRIuso5d7ikEyV9M/p3lfubG3ZKOs05d4ykcZKmmdmJkn4mabZz7kBJWyRdMnBTRJqukLQ85jH3NrdMds6Ni1kfOKs+mwnBiZ0g6UPn3ErnXLOkeZLOGuA5YQ85516StLnL4bMk/T768+8lnd2fc0JmOOfWOufeiv7cqMi/TEeJ+5sTXMTW6MP86C8n6TRJD0ePc38HKTMbLWmGpN9GH5u4t7kuqz6bCcGJjZIUjHkcih5D7tjLObc2+vM6SXsN5GSQPjMbK2m8pNfF/c0Z0f9dvlTSBknPSPpIUr1zriU6hM/nwWuOpGsktUUfV4h7m0ucpIVmtsTMZkWPZdVnc95AXhzIBs45Z2asFTiImVmJpD9LutI5F44UlCK4v4Obc65V0jgzGy7pEUmHDuyMkAlm9llJG5xzS8xs0gBPB33jU865WjOrlvSMmb0X+2Q2fDZTCU6sVtI+MY9HR48hd6w3s5GSFP19wwDPB3vIzPIVCcAPOOf+Ej3M/c0xzrl6Sc9LOknScDNrL+Lw+Tw4nSzpc2a2SpGWw9Mk3Snubc5wztVGf9+gyH/AnqAs+2wmBCf2pqSDot9SLZB0vqTHBnhOyKzHJF0Y/flCSfMHcC7YQ9Eewt9JWu6c+2XMU9zfHGBmVdEKsMxsiKTTFen7fl7SudFh3N9ByDl3nXNutHNurCL/jn3OOfcVcW9zgpkVm1mg/WdJZ0hapiz7bGbHuCTM7DOK9Cv5Jd3rnLttYGeEPWVmD0qaJKlS0npJP5D0qKSHJI2RtFrSec65rl+eQ5Yzs09J+pukd9TZV3i9In3B3N9BzsyOVuTLM35FijYPOeduMbP9Fakelkv6h6SvOud2DtxMkY5oO8TVzrnPcm9zQ/Q+PhJ9mCfpf5xzt5lZhbLos5kQDAAAAM+hHQIAAACeQwgGAACA5xCCAQAA4DmEYAAAAHgOIRgAAACeQwgGAACA5xCCAQAA4Dn/H+MMxBK0tGqOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mod = layered_earth_model(v, z, theta, td)\n",
    "plt.subplots(figsize = (12, 4))\n",
    "mod.plot_model(model_span)\n",
    "mod.plot_layout(sources, receivers)\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def model(params, config):\n",
    "    v = params[0:2]\n",
    "    z = [config[2], params[2]]\n",
    "    theta = [config[3], params[3]]\n",
    "    sources = config[0]\n",
    "    td = params[4:len(sources)+4]\n",
    "    receivers = config[1]\n",
    "    ic = np.arcsin(v[0]/v[1])\n",
    "    x_span = [min(min(sources), min(receivers)), max(max(sources), max(receivers))]\n",
    "    times = np.array([np.zeros(len(receivers)) for i in range(len(sources))])\n",
    "    for j in range(len(sources)):\n",
    "        h = (z[0]-z[1] - sources[j] * np.tan(theta[1]-theta[0])) * np.cos(theta[1]-theta[0]) # Thickness of layer one at the source point\n",
    "        t_direct = [abs(receivers[i]-sources[j])/v[0] for i in range(len(receivers))]\n",
    "        t_refracted = [(abs(receivers[i]-sources[j])/v[0])*np.sin(ic-np.sign(receivers[i]-sources[j])*(theta[1]-theta[0])) + (2 * h * np.cos(ic))/v[0] for i in range(len(receivers))]\n",
    "        stacked = np.dstack((t_direct, t_refracted))\n",
    "        # print('j: %d, td[j]: %.3f' %(j, td[j]))\n",
    "        times[j] = stacked.min(2).flatten() + td[j]\n",
    "        # print(times[j])\n",
    "    return times.flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fun(params, config, obs_times):\n",
    "    return model(params, config) - obs_times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "config = [sources, receivers, 0., 0.]  # x locations of sources, x locations of sources of receivers, ground elevation at x = 0, ground slope\n",
    "starting_params = [750., 1800., -7., 0.09, 0., 0., 0.]\n",
    "picking_noise = 0.003\n",
    "theoretical_times = mod.two_layers_forward_refraction(sources, receivers)\n",
    "if picked_times == '':\n",
    "    obs_times = theoretical_times.flatten() - (0.5 * picking_noise * np.random.laplace(size=len(theoretical_times.flatten()))) + 0.5 * picking_noise * np.random.laplace(size=len(theoretical_times.flatten()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-2-13a568e85554>:83: UserWarning: marker is redundantly defined by the 'marker' keyword argument and the fmt string \".k\" (-> marker='.'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: marker is redundantly defined by the 'marker' keyword argument and the fmt string \".k\" (-> marker='.'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: marker is redundantly defined by the 'marker' keyword argument and the fmt string \".k\" (-> marker='.'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: marker is redundantly defined by the 'marker' keyword argument and the fmt string \".k\" (-> marker='.'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string \".k\" (-> color='k'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: marker is redundantly defined by the 'marker' keyword argument and the fmt string \".k\" (-> marker='.'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string \".k\" (-> color='k'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: marker is redundantly defined by the 'marker' keyword argument and the fmt string \".k\" (-> marker='.'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string \".k\" (-> color='k'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7faa31f6b2b0>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAD4CAYAAAAejHvMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAA4Y0lEQVR4nO3dfZQc1X3n//dFsiRwjGyY8R4WQSTW2EYyIGwJo4A1D0g2+AHWGMs4EAkfJwRhcmxngwIko+kZ2fGitZb4HGMIwUayE/+MrPzsJdkNgngGgQnGM9giSOinBAk5IMhqRjz87MWOEfruH1U9qmnNTNdDV3V11ed1zpzprq6euj1VXfWte7/3XmdmiIiIiIjIEcc0uwAiIiIiInmjIFlEREREpIaCZBERERGRGgqSRURERERqKEgWEREREakxvdkFqNXW1mZz585tdjFEREREpOAef/zxUTNrn+i13AXJc+fOZXh4uNnFEBEREZGCc879bLLXlG4hIiIiIlJDQbKIiIiISA0FySIiIiIiNRQki4iIiIjUUJAsIiIiIlJDQbKIiIiIjLN+PQwOjl82OOgtr6dSqaRSpqwpSM6hJAemiIiISFKLF8OKFUfikcFB7/nixfXf29fXl27hMqIgOYeSHJgiIiIiSXV1webNXvyxdq33e/Nmb3lZKEjOoTwemKrdFhERKZeuLli9Gtat835PFYdUKhWcczjnAMYet3LqhYLknIpyYGahrLXbujkQEZGyGhyE228H6Of224++HgZVKhXMDDMDGHusILmk0tzx1QOzp4e6B2YW8li7nYWy3hyIiEg2sq6MCRu7VK93mzcD9I7FAM2ORzJVjfTz8vOe97zHWoX374umt7e37joDA2Ztbd7viZ43U0+PGXi/y6L6/+/pyc9+EBGRYohzzb/llqNfHxjwltcTNnZZtux+g04DAj+dtmzZ/XXfGybWyQtg2CaJSZ351eJ5sWjRIhseHm52MUJxzhH1/xfmPevXezWVXV3eHV+lUmFwEIaGYM2aJCVOpnpXuXq1V7tdhprkqrVrvdSXnh7o7292aUREpEiiXl+DtbzbtlXo6KiEbuFNK3ZpVc65x81s0YQvThY9N+sn7zXJvb29NXdV3k/YuyYi1j5HXT8tea7dTptqkkVEJG1RW2qr1yboq3ttyjp2aSVMUZMcKnAFLgJ2A08DN07w+kzgHv/1x4C5/vIrge2Bn8PAwqm2lfcgOSjsQZPk4Ix6YKbVxJGkaSdP24iqzDcHIiKSjbiVMXFSIOMEvK2UPhFVoiAZmAbsAU4DZgBPAPNr1rkOuMN/fAVwzwR/50xgT73tFTFIjvqeLIPq6vbyII8BaR4D9ziK8jlEpHXovBNO8FrX29sb6trnxQmdBgcM+vzfnZNez4P7ohonaF94kgbJS4Ctgec3ATfVrLMVWOI/ng6MgpfvHFjnz4Av1tteKwXJcYLLtNMt0grcg8J+7jgnSKU2pCOPNyAiUmw674QTJ4AN/i+Buv/bOIF4WSQNki8H7go8/x3gqzXr7ADmBJ7vAdpq1tkDvGuSbVwDDAPDp556aib/lGaJGlinXfMcdhtx1o/7pSzjCBpZ0A2IiGSt7OedNK75ZskC67Lui8k0PUgG3gs8WW9b1mI1yVlI6wuWVUpHlI4FwfX1JU6HbkBEJGtlPu9kUdEVJU4o876YTNPTLYBbgZvrbcsUJCeWxzzpsF9KNQelSzcgIpK1op130qq4irt+FEXbF42SNEieDuwF5nGk496CmnU+w/iOe5sDrx0D7AdOq7cta4EgOY8dEYJlqn6Bo5QprS9xFh0Lkn72slBuoIhkrYjnnbQrldIKkou4LxolUZDsvZ8PAv/sp1H8ib+sH7jEfzwL+C7eEHA/DgbEQCfwozDbsRYIkvN4oCUtU1p3xlE7FmS9jbyM6pGFPN7ciUix5X3o0Dx2vtcwrtlLHCRn+ZP3INksn00WWZYpzugWYWqG4959R817rsrLySuOVj0ZiYg0UpJKorz14Wl1eaxEDENBcgrymPwepUxZB1npd0DsMzD/dzonrzyd7Fr1ZCQi0mhxK4nS6sMTlKfKlSzksRKxHgXJDZbHgyBqmfIeZKU1gkaRJmrJ43EoItIMYSuJshgytewtfXmsRJyKguQGymNwGbdMeQ6ywp6w0s57znoM6qhBdaudjEREGi3tlLuoHcTzGCdkJc9xxWQUJDdQHu8Qk5Sp1YOsqCevJFNzptU0F7dM1ZPR0qUPtszJSEQkjLxUlNRuY6LnU72nlYLFpFr15kBBskyo7F/iqGMxp9XJI85JXqN6iEiRhT3fLlt2vz/UaPB822nLlt1f971RzoVxrpetXgkVVR4rEcNQkCxHadU7vkaIe3OQ1nBBUceTNsu+RlxEJIm0J+GI+54oogS9ZayEalUKkuUorXrH1yhp3uHHDWCjlKlIHRBFpPjCnHeS1AqH3UZcUYLeMldCtSIFySIBad/hx0npiNvxxCxKbXV2Pbo1+6GIBIU5hyRJIzNL72Y+atBb9kqoVqMgWcSX1R1+3FqHqHnSZtFrT9JqxlSetEjryGKmujg350kqDMKI87kV9BabgmQRX5Ynu7DpE0nLlP5ELelf4JQC0li6qKenKP/ba64xmz17fIXB7Nne8nrSzhdOMx1OqRBSS0GySMby3Gkj3Q6I4Wc/TNL5MGyZyipOIBA1+CtKsBhVUYKsgQGz44/3AuOlSx+02bO9540c6SfOe7I4d+b5/CzZU5DcZGW9mJRVUS6iQVEvcGFqkleuvHtsJI9qBx04YCtX3j3pe5IE1mWreY4aCEQ9bot4nIdVlCBrYMDs2GO9SODYY6f+HEn7NYRZL8tjqmzDs8nkFCQ3WZkvJmVUxJuiqBe4qOM9h03PiBNYV6WV0pHn/R01EIgbWLd6sBhHnoOsKBNxHHec9zmOOy7d9Kgwsvoulfm4laMpSM4BfSml6OKObhE3kMtiCtowtdVJO16mJe45J+r+yHOwmIZbbjHbsGH8/3bDhubfFEU9br0bzpf8n76xx2ndcOaFKq2kloLknCjbxUSknuoFKurU2mG/S0k6H6ZZIx4sXxriBgKqSa5vwwYz57zfEz1vtLSmZ77mmiM5yNX1jz8+XMe9Vk5dynPLjzSHguQcKOPFRGQqSQO5qIF1tFE9os2AaBbvJjhPzdbKSQ4n65rktI7bpB1nRYpCQXKTlfViIjKVLAK5oKgBabTa6uhBdZwyVbeXhiKMbpHHIR6TSuu4DWrlmmGRpBIHycBFwG7gaeDGCV6fCdzjv/4YMDfw2lnAo8BO4Elg1lTbKmKQnMeLiUgrymIChOrfDJs+EbWZO+upd8sUAGVVIZGkZTDM/og7kkTctB+RMksUJAPTgD3AacAM4Algfs061wF3+I+vAO7xH08H/gk4239+IjBtqu0VMUgWaUVlvbmL2hEvSUe/OFPvZjHDYisH1mmntiUNxNPaf3ntQCqSd0mD5CXA1sDzm4CbatbZCiyxI4HxKOCADwJ/VW8bwR8FySL5UNY0oSxqq6PW+CWpfU57djSz/AXVaaZCxB3FpSqtIDntm9qy3jRL8SUNki8H7go8/x3gqzXr7ADmBJ7vAdqAzwHf8oPonwBrJtnGNcAwMHzqqadm9X8RkTpavcNpni/sUQK5qLXPSSd+yKK2Oi1Rj9kkx0jYz53kJicvNyBlvWmW4mtmkPxHwDP+4+P83OQLp9qeapJF8qWVhy7M64U9zs1H3HzTsIFcksA6LykdwWHNzCzUsGZJ0hTipELESbHJi1a/aY4rzzfbklwz0y2uADYF1usBbphqewqSRfKjCBfFvH2GJIF7VsPMhXlP1rXVYf7uNdeYzZ49/n87e3b9sX+j3IBk0akuz0FZK980x5XXm+205fk4bKSkQfJ0YC8wjyMd9xbUrPMZxnfc2+w/foufZnGc/3f+AfjQVNtTkCySD0W6MOTpwh73whM32I9TY5tFukVa26j+n9KaoCZOmaJuI6/fvbzdcGapjJ89r8dhoyUKkr3380Hgn/00ij/xl/UDl/iPZwHfxRsC7sfAaYH3XoU3/NsOYH29bSlIFsmHotQiFOHiFvVilXTfRQ2s85bSESXgzWqc66jHYd6O27IETFPJ0812VqIch616zUgcJGf5oyBZRBqlKBf2qBefrD93WrXVSaYVD5u7nSRfOM6U0RM9n0yegrJWDYAaJW83LVmK0wLSSsMQKkgWkVIq84U9jxf1JFMhhwmqV668e6wmuDqCBBywlSvvTqVMYYPkJLNL5mn/lVVRbrbjiNsC0koT2kwVJB+DiEhBrVkDXV3jl3V1ecuLrqsLVq+Gdeu837X/h2bYswc++lEYHITe3l4GB73ne/Y05u8///zJwArgQX/Jg8AKf3nj9fX1RX5PpVKpu87gIKxYAZs3Q3+/93vFCm+5ZG9oyNsH1e9QV5f3fGioueWKYv36o4+fwUFv+WTiHIfV8w6szc15J5HJoudm/agmWUQkuTzWRA4MeEOyzZ7tdaybPXv8kG1TIaWOfnFqn6NuI2pKR9IJS0RqxakNj3ocJsnvbyaUbiEiUh55bh4eGDA79ljv6nPssekMt2YWLajO2xBwcT5HmVOLJJwkN85hjsMk+f3NNFWQrHQLEZGCybp5OEwKQZBz439P9XerFys4UqkTZnu9vb2hyxOliThOmSqVCt3djtHRfmAto6P9dHe7Sd8Tp2l88eIjTeGVSmWsqXzx4snfI+WSdgpWEdJSjjJZ9NysH9Uki4hkL4vpmb3Uhpf8n76xx41MbYgj7tjKUcoUpSY5bgpIK3aakuxEPT6StOTkPcUiCKVbiIjIVJKkaIQNFoPTRuM3x9abNroqrYtuks8dZwg4QjZDxw148zRsnGQjzHGYNBUizZvUZlOQLCIidUWpUY1Ty5RkuLW0ZJHLG7cjXlaTokhrCxPALlt2f6BlovrTacuW3d+wbbQqBckiIhJKFtMzm7VWc2wzRE0BadVOUzJeXqeRL/L3daogWR33REQKLEqnusFBuP12gH5uv33q8VDjdC6LW66yCY5Pu21bR6jxaQvZaaqEwo69XalUcM7h/N6v1cdpfa/K+n1VkCwiUmBhL7qrVm2ku3uE0dEuoJfR0S66u0dYtWrjhOsHR1OoTgzS6NEUkgbirSpOwBucOKc6skdZJs4po6xGfik7BckiIi0krRqdBQuuZmCgHTMvKjUbZGCgnQULrp5w/WrgtmIFHD5cGav5bOSwUsFAHCjNsGZJZ4osa61fXgRv7qr7Yqqbuzi1wkluIHV8RDBZHkazfpSTLCIyOULkEyYZuinsNqrSHk0hjzMHFk2R802bIUl+eNjvXvBv9vb2Kgc9AdRxr3E0q5GINFOUADbO+mbRhzZLO4DVsGbpinOMlFXU70ZaMywGt6EbyGSmCpKVbhFRWZv/RKTxwjZ75rGTTrBzWX8/oTqXxVHtTNjTQ93OhCJpC5Pjn2SGxWq+cJjUibRn0BNUkxyH7t5EmqNoLTnEqMGL+p60mtKz2BdJJvqQqSVNySmCNIdbi1KTHPc4VyzSGCjdovHU/CeSvbwHTVEvulkEya2saDdFeVWmYyoo7OeOekORZIbFsAFv3s+FrURBcoPp7k2kefL8/Qtz0U1ag1emmj7JhoLkxr4nixkWdQPZOImDZOAiYDfwNHDjBK/PBO7xX38MmOsvnwv8Etju/9xRb1t5D5J19ybSfFm15KRdM1zW4ETypSg3XmE+R5Yjv0SR55v/oksUJAPTgD3AacAM4Algfs0611UDYOAK4B47EiTvqLeN4E/eg2TdvYk0V5YXk7RrhhUkS6tohWtfFjepadxQqPKtuaYKksOMbnEu8LSZ7TWzXwPfAS6tWedSYJP/eAtwoat2wy6YpIO8i0h8SUZUCDsSRNRB+qPOfBWnJ3tZZ54r6+fOo6xHdkpr5JY8TqeuKcVzbLLoOXDSvxy4K/D8d4Cv1qyzA5gTeL4HaMOrSf4/wE+BbcD7JtnGNcAwMHzqqaemf9tQMK1why/SCEmOdULWGq1cebfBAYNOv1a40+CArVx5d0O2EafWqKw1TWX93HmVditO8Ptd/S7V+34n6VSnSTjELHm6RZIgeSZwor/sPcCzwPFTbS/v6RZ5pAuJSH1hg2SzaMM3xemkEyfYKGvOYlk/d17F7Q8QJk0hyUx1ZuG/4zqmJChpkLwE2Bp4fhNwU806W4El/uPpwCjgJvhbDwKLptqeguR49KWXskm7k07YYCDuTWqcYKOsQ0+W9XPnTdRZ5KLWDHvf106/JadvrEWnFadTl9aRNEieDuwF5nGk496CmnU+w/iOe5v9x+3ANP/xacB+4ISptqcgOT596aVMolwQo65fDQaWLn0wlUH9VZMcXlk/d94kGfs3as1wmrXVwXLpmBKzhEGy934+CPyzn0bxJ/6yfuAS//Es4Lt4Q8D9GDjNX/4xYCfe8G8/AT5Sb1sKkuPRl15aWZozX0VdP+2aYeUkh1fWz51Hy5bdH8jTr/502rJl90/6njg1w1Frq6PSMSW1EgfJWf4oSI5OX3ppdWED2CTpE2ED8TidA6PcpMb5+2XtnFvWz513UW9Q49xAptWpTseU1FKQnJK8DMCuL73kTZmmZ9ZNqrS6NPP743aCDb5fN5CSJgXJKYlzkc5LYC1ST9rDrWUx81UWF1BdpKXVZZG6lFbNsG5SJSkFySnJovZLQbU0QpIUgjgXuDQ71VVFHVJqouciEv2GM8xIFWZmF19stmHD+GUbNnjLG0l9ciQJBckNlPW873lpUi6KstYsxg0WozSV5nV6Zl1ApQyinneidsSLc9Oc5U2qRneSuBQkpySLzkaqeW6sLE7aea29jBssxrn45O241QVUii7qeSfO8GxxziFZ3KTqRliSUJCckrTSLbKukStbYF3mk3aUYDHq8E1xmmOzkNd9IdJoUY/1OMOt5W0SnLxWSkjrUJCckjyN7Rp3/bDviTP1bp5lUbOYt9rLqBfEqDVNWXTSiUoXUCmbqOedKOsnqUkOOzFPVHlMb5PWoiA5R9IYGiuLPOk4TXO1Zawnq5Nd0WqSo3ZgC7v/4tQM563WVhdQKZO4Nclh1o9zw5n0uiGSBQXJLSyt8WbjBNZxmuailKso+cJZ116G+d/GmS0rKMpxmLcadJEyiNqSE/U8lXQM4zylYIkEKUgukfTzpPsMzP/d+FE9itA0l2QbeUzhiSJvNckiZRE1IM3iXJi0lVMkCwqSSyTNICtqTXLcE2Tcmsg0TrxZN9fnccSUsLKomRKR+vI4dGgeyyRipiC55WQdPKSV0xqUViAeZxtVeZyMopUnqIl63KpTnUjj5L3WVkGy5JWC5BaTRfAQNaCJM7pF1Oa/lSvvHhtmrJozCwds5cq7G7aNoCxSQMJcoJLmC7fyxUfpGSKNl8dzQl6CdZFaCpJbUNrBQ9Yd2MJ0JIkT8MYJrKuiXEjipoCkNXJIkYbkU0c/kcbKY5AsklcKkltU2sFDnodCi3KSj5KiEbVJ8pZbzDZsGP/3N2wIH4ymlWZSlFQF1SSLNJ5qbUXCmypIdt7r+bFo0SIbHh5udjGabnAQVqyA1avh9tth82bo6mr8dtauhXXroKcH+vsb//fjbqNSqVCpVFLdhnOOesf/Bz5wH/ff/37gvwB/DnwO2MD7338/W7deNOF7li9/gH/4hz8DHgws7WTZspt54IHlDfsMWR0jaamWv1ru2uciUi7r18PixeO//4ODMDQEa9Y0r1xJvfbaazz33HP86le/anZRSm3WrFnMmTOHN7zhDeOWO+ceN7NFE75psui5WT+qSc6uljDPNclhBWt6qznDYWt6CTnTYNSa5KjpE8FtVP9PYT9DK6cqaHQLkYmV9bsRNUWvVezdu9dGRkbs8OHDzS5KaR0+fNhGRkZs7969R71G0nQL4CJgN/A0cOMEr88E7vFffwyYW/P6qcAvgD+qty0FydmcILPOSU5rGxs2mDnn/Z7o+VTCBMlxxof23tPp50r3jeVMT/aeuJ9BqQoixVSUdKo4koxwlFdPPfWUAuQcOHz4sD311FNHLU8UJAPTgD3AacAM4Algfs061wF3+I+vAO6peX0L8F0FyfmR90k1omwjbi1smCDZLP5JO2wtb5zPUOaLqEgZlPkmuJVbyCYyUWAmzZFGkLwE2Bp4fhNwU806W4El/uPpwCiM5Tv/Z+C/ARUFyZKWsCfVqB334o4PHSewjnJhKGtzrEiZFC1YrCdqK1yraHaQ/NJLL9ltt91mZmaDg4P2oQ99KNPt33333bZ///6x55/+9Kdt586dkf9OI8oeNUg+pl6iM3Ay8Gzg+XP+sgnXMbNDwCvAic653wD+GOibagPOuWucc8POueGRkZEQRRI5YnDQ67jW0+P9HhycfN1KpRK8uRt7PFknwaGhIx3Jent76eryng8NTV2eage03t7DbN7sPZ+qXNXPAP11PwN4nVhqO7d1dbV25xaRolq//ujv9OCgt3wyUc8JRdDRUaGtbZCBgXagl4GBdtraBunoqDS7aE0RpfP6VF5++WW+9rWvNeRvTebQoUOTvrZx40aef/75sed33XUX8+fPT7U8jRImSE6iAtxqZr+YaiUzu9PMFpnZovb29pSLJEUSDEiPOaYSKiCNIhiMVk9Y9YLRYGBdqVTqBtbBzwC9Df8MItJcixeP/05Xv/OLF0+8flnPCcFzJxCqUqLI+vqmrF8M7cYbb2TPnj0sXLiQG264gV/84hdcfvnlvPOd7+TKK68cqzR6/PHH6ejo4D3veQ8f+MAHeOGFFwDYvn075513HmeddRYf/ehHeemllwDo7Ozkc5/7HIsWLeIrX/nKhO/fsmULw8PDXHnllSxcuJBf/vKXdHZ2Uh3F7L777uPd7343Z599NhdeeCEAP/7xj1myZAnnnHMOv/Vbv8Xu3bsb8n+IZbIq5kBtW+x0C+BhYJ//8zLwInD9VNtTuoVEkWTGvbw04SWdcU9E8i9KjrHOCfk5PzdC3HQLQvabqeeZZ56xBQsWmJmXsnD88cfbs88+a6+//rqdd9559vDDD9uvf/1rW7JkiR04cMDMzL7zne/Ypz71KTMzO/PMM+3BBx80M7Oenh777Gc/a2ZmHR0dtnr1ajOzKd/f0dFhQ0NDY+WpPj9w4IDNmTNnbMSJgwcPmpnZK6+8Yq+99pqZmT3wwAN22WWXjZU963SL6SHi6CHgdOfcPGA/Xse8365Z515gFfAocDkw4G/4fdUVnHMV4Bdm9tUwwbtIGBPV6HZ1TT7ObnAczmrNcLPH4fTGTvbGTw4zdrOItJ6uLm9M8+pY6FONBa5zQuNSDVpNpVIZV4PsnAO8dL9G/U/OPfdc5syZA8DChQvZt28fb37zm9mxYwfLl3vH3euvv85JJ53EK6+8wssvv0xHRwcAq1at4uMf//jY3/rEJz4BwO7duyd8/1R+9KMfsXTpUubNmwfACSecAMArr7zCqlWr+Jd/+Recc7z22msN+dxx1E23MC/H+Hq82uJdwGYz2+mc63fOXeKv9nW8HOSngT8EbkyrwCJBlUoF59zYiaT6eLKTSdRmTxGRoDj5xdV1wvadkPKK2m8mjpkzZ449njZtGocOHcLMWLBgAdu3b2f79u08+eST3H///XX/1hvf+MaxcsZ5/0R6enro6upix44d/O3f/m1TJ2EJlZNsZv/LzN5uZv/JzL7oL1trZvf6j39lZh83s7eZ2blmtneCv1Exsy83tvhSdlFPKNUctxUrvBnu8jbDW29vb7OLICJTiHOjHcwx7u8nUo6xzgmS1Jve9CZ+/vOfT7nOO97xDkZGRnj00UcBb5bAnTt3Mnv2bN7ylrfw8MMPA/Ctb31rrFY5zPun2v55553HQw89xDPPPAPAiy++CHg1ySef7I0PsXHjxhifuHHS7rgnkjvBZs/Vq/MTIEN5mxhFWkWcG+0kHdJ0TgivaP+rRt0gnXjiiZx//vm8613v4oYbbphwnRkzZrBlyxb++I//mLPPPpuFCxfyj//4jwBs2rSJG264gbPOOovt27ezdu3aSO+/+uqrufbaa8c67lW1t7dz5513ctlll3H22WePpW6sWbOGm266iXPOOWfKUTOy4PKW67Ro0SKr9noUiaJSqYQ6SVZrdVav9po981STLCKtYe3aI/nF/f3NLk15BfuZVPO3m93PpNauXbs444wzml0MYeJ94Zx73MwWTbS+apIlN+Lm+lVFCZDjNHvmSdL/lYjEp/zi/FA/E0mTgmTJjSxOdkUZhzP4v6pUKrowiGQk6Y120VICmm3btgqjo110d48AfXR3jzA62sW2bZVmF00KQOkWGQg2B1XlrTkoL5QKEV71fzU62k9b21r9r0QykPR8XtYh3dKW5/QXpVvkh9ItckjNQeHluVNd3lT/V7BW/yuRjGhK+Pwp4xTekg0FyRnI+7BjeVI92S1duk0nuyl440N3sW7dCNDPunUjONelplyRHIo6nruEF0x/6e09XDf9Rf05JAoFyRlRDWl9wZPdQw91tmynuix0dFRoaxtkYKAd6GVgoJ22tkE6OirNLpqI1MhigoiyCvYzqVQqdfuZqGVXolCQnBH1hq6vKJ3qsqD/lYhI9PQXteyO97u/+7s89dRTk75eqVT48pfjzQO3ceNGrr/++qOWP/jgg2NjKAPccccdfPOb34y1jbQpSM5AUYYdS9urr1bo7h7fJNnd7Xj11UpzC5ZDwQtDdcB55UWK5J9m0Gu+PLfsZp0OctdddzF//vx0/vgkaoPka6+9lpUrV2ZahrAUJGdAtX7hqEkyHv1/RFqHvq/Nl+eW3TTSQfbt28c73/lOrrzySs444wwuv/xyXn31VQA6Ozupjih233338e53v5uzzz6bCy+88Ki/85d/+ZdcfPHF/PKXv+Sv/uqvOPfcc1m4cCG///u/z+uvvw7A3Xffzdvf/nbOPfdcHnnkkQnLcscdd3DrrbeycOFCHn744XG11Z2dnXz+859n0aJFnHHGGQwNDXHZZZdx+umn86d/+qdjf2ei7b/++utcffXVvOtd7+LMM8/k1ltvjf9P8ylIzoB6Q4uIiDRfsGX3mGMquWvZTSsdZPfu3Vx33XXs2rWL448/nq997WvjXh8ZGeH3fu/3+Ju/+RueeOIJvvvd7457/atf/Sp/93d/x/e//3327dvHPffcwyOPPML27duZNm0af/3Xf80LL7xAb28vjzzyCD/84Q8nTOOYO3cu1157LZ///OfZvn0773vf+45aZ8aMGQwPD3Pttddy6aWXctttt7Fjxw42btzIwYMH2bVr14Tb3759O/v372fHjh08+eSTfOpTn0r2T0NBsuSUmiRFRKTRgi27fX19uWzZTSMd5JRTTuH8888H4KqrruKHP/zhuNd/9KMfsXTpUubNmwfACSecMPbaN7/5Tf7+7/+eLVu2MHPmTH7wgx/w+OOPs3jxYhYuXMgPfvAD9u7dy2OPPUZnZyft7e3MmDGDT3ziE7HKeskllwBw5plnsmDBAk466SRmzpzJaaedxrPPPjvp9k877TT27t3LH/zBH3Dfffdx/PHHx9p+0PTEf0EkBWqSbCxNaCMiMvH5rqsrX3nJtekgjShfta/PZM+ncuaZZ7J9+3aee+455s2bh5mxatUqvvSlL41b7/vf/36yQvpmzpwJwDHHHDP2uPr80KFDk24f4IknnmDr1q3ccccdbN68mW984xuJyqKaZJES0DTWIo2nm/nWk/cxq9Pq6P+v//qvPProowB8+9vf5oILLhj3+nnnncdDDz3EM888A8CLL7449to555zDX/zFX3DJJZfw/PPPc+GFF7JlyxYOHDgwtu7PfvYz3vve97Jt2zYOHjzIa6+9dlTKRtWb3vQmfv7zn8f+LJNtf3R0lMOHD/Oxj32ML3zhC/zkJz+JvY0qBckiJRDMc+vrO6b0wx6JNEJfX1+ziyAR5b2DeFod/d/xjndw2223ccYZZ/DSSy+x2puudUx7ezt33nknl112GWefffZRqRIXXHABX/7yl/nQhz7EW9/6Vr7whS/w/ve/n7POOovly5fzwgsvcNJJJ1GpVFiyZAnnn3/+pFNxf+QjH+F73/veWMe9qObPnz/h9vfv309nZycLFy7kqquumrCmOSqXtznkFy1aZNWelpIfSZvrK5VKbk5CZbZ2rZfn1tPj1VKISHzOOfJ2DZXwstp/u3btmjRgzMK+ffv48Ic/zI4dO5pWhryYaF845x43s0UTrR+qJtk5d5Fzbrdz7mnn3I0TvD7TOXeP//pjzrm5/vJznXPb/Z8nnHMfjf6RJA+SDkujGpfm0jTWIo2R9+Z6CU8dxKWeujXJzrlpwD8Dy4HngCHgk2b2VGCd64CzzOxa59wVwEfN7BPOueOAX5vZIefcScATwH80s0OTbU81yflVDYxXr/Y6E9Rrrg/WPlfv2NVZrDmCeW7d3Y6BAVPKhUhCqkmWMJpdkyxHpFGTfC7wtJntNbNfA98BLq1Z51Jgk/94C3Chc86Z2auBgHgWoLNJC4s6LM3OnRvp7vZqLAGc66K7e4SdOzemX1gZRxPaiIiIRBMmSD4ZeDbw/Dl/2YTr+EHxK8CJAM659zrndgJPAtdOVIvsnLvGOTfsnBseGRmJ/ikkE1FnKdq06WoGBtppaxsE+mhrG2RgoJ1Nm67OorgSoGmsRRpPzfUSllocmi/OPkh9dAsze8zMFgCLgZucc7MmWOdOM1tkZova29vTLpLEEHdYmmrtM6xt2KDokoxyJ0UaQ9+l8omzz2fNmsXBgwcVKDeRmXHw4EFmzToqBJ1SmMlE9gOnBJ7P8ZdNtM5zzrnpwGzgYE0BdznnfgG8C1DScYuZqrl+qsC3Wvu8dOk2br+9I3eDtouIiEwl2L+mr69vbKz5sP1r5syZw3PPPYdayptr1qxZzJkzJ9J7wgTJQ8Dpzrl5eMHwFcBv16xzL7AKeBS4HBgwM/Pf86zfce83gXcC+yKVUHIhzixFwdrnrq6OmufplVVERKRRqqM7bd7sPQ9ey8J4wxveMDbds7SWuukWfg7x9cBWYBew2cx2Ouf6nXOX+Kt9HTjROfc08IdAdZi4C4AnnHPbge8B15nZaIM/g+SUOosVh5qVpQx0nMtEtm2rMDrqdTyHPrq7Rxgd7WLbtkqziyYp02QiIlKXhrqSotJQlRKWJmMqpsSTiYgEqbZFRIoi6URJUg7V/jXQH2p0JykGBckSmWbPKwfNLCZloKZ0qSeYg9zbezj06E7S+pRuIZGp6b18tM+l6NSULpMJpuRUKSWnOJRuIYmpVlFEikpN6TKV4GRMVZqMqRzCDAEnQqVSGQuIVatYPppZTIoq2JS+bdthOjo0VKWIeJRuIZEpSBaRolBTuki5Kd1CGkq1iuWwfv3Rzc6Dg95ykaJQU7qkSSmJrU1BckklCYD0pS8HDY0lIpKMRoNqbQqSS0oBkNRTnSFxxQqv53+UPE3dSEmz6RgUkaQUJJdUkgBIyqOrC1av9obGWr06/PGh2hNphmALWfUYVIqQZE2jQRWHguQSixsASXlUh8bq6UFDY0nuqYVM8qBSqWBmYx3cq48VJLceBcklVg2Ali7dpgBIjhIcGqu/n7qzTKn2RJpNs+eJSCNpCLiSCgZA3d2OgQFTyoWMk2RoLA0TKM2k2fOkmYLnzuocAxpWML+mGgJOQXJJBb/E1YBGX2JpFAXJ0izVCoDR0X7a2tbqxl8yF6yE6uo6+rnki8ZJlqO8+mqF7u7xTePd3Y5XX600t2BSCBpLW5ohGIz09h6umyIkkgZ1jC8O1SSLav1EJBNRU3jSXl8kTUr7aQ2qSRYRkaaLOvpE1PU1e57kRdyRgdTROV8UJIuaxkUkE1GbodVsLa0o6shAQRpjPl9CBcnOuYucc7udc087526c4PWZzrl7/Ncfc87N9Zcvd8497px70v/d3eDySwPozlWaIcnU6NK6oo7PrvHcpdUMDY2/mave7A0NNbdcEl3dINk5Nw24DbgYmA980jk3v2a1TwMvmdnbgFuBW/zlo8BHzOxMYBXwrUYVXERaW7ApvTpEkiZ+KL6ozdAaz11aTdS0n+XLH8C5rpox5rtYvvyBlEsq9YSpST4XeNrM9prZr4HvAJfWrHMpsMl/vAW40DnnzOynZva8v3wncKxzbmYjCi4irS3YlN7Xd4ya0ksgajN0cP2HHurUaBVSSDffvJy2tkEGBrwO9AMDRlvbIDffvLzJJZMwQfLJwLOB58/5yyZcx8wOAa8AJ9as8zHgJ2b277UbcM5d45wbds4Nj4yMhC27TECpE9JKqk3psFZN6SUQtRlazdZSBsEKA+hThUGOZNJxzzm3AC8F4/cnet3M7jSzRWa2qL29PYsiFZaS/qVVeNNYd7Fu3QjQz7p1IzjXpRu9AovaDK3x3KUsVGGQT2GC5P3AKYHnc/xlE67jnJsOzAYO+s/nAN8DVprZnqQFlqOpA5TkQdTjsKOj4jcxtgO9DAy009Y2SEdHJe2iSouoVCqY2dg47tXHupGSook7ZJykK0yQPASc7pyb55ybAVwB3Fuzzr14HfMALgcGzMycc28G/idwo5k90qAyF17UYGPnzo10d3u1cADOddHdPcLOnRvTLahIQNQxbdWULiKSbMg4SVfdINnPMb4e2ArsAjab2U7nXL9z7hJ/ta8DJzrnngb+EKgOE3c98DZgrXNuu//z1oZ/ioKJGmxs2nT1WC0c9I3Vzm3adHVWRRaJPKZtsOm9Ola3Jn6QyWg8dykqVRjkl6alzqlqYLx6tdf0EiaJX1NgSh7oOBQRyU6lUlEKUgKalroFRR1AX2OJSh4or05EJFvqsJ8eBck5FSXYCOYzbdvWoXwmaQrl1YmISJEoSM6hqMGG8pkkD7I+DtW82Nq0/0Ti84bQHD88onNO36sGU05yDq1f73XSC6ZYDA56wYY6NYl4nHPk7fwl4Wn/iTSGvkvJKCe5xQR7/VfvCtXrX0RERCQ7CpJzTgn5IkeoibG1af+JNJ6GR0yP0i1yTs0oIhPTd6O1af+JSB4o3aLFqLZFREREpLmmN7sAcrTgwOCqbRGZmJoYW5v2n4jkndItck5BsoiIiEg6lG7RwlTbInLE+vVHjxc+OOgtb8T6IiJFp9TN8BQk55wOZpEjFi8eP7FOdeKdxYsbs76ISNFp1KzwFCSLSMuozuK3YgWsXXtkZsrgxDuTrd/Rsa3u+hKPbuZFpIgUJItIS+nqgtWrYd0673e9gLe6/kMPdYRaX6JTzZRItqKmkmnUrHgUJItISxkchNtvh54e73fthWKy9aE/1PoiInkXNZWsUqlgZmMDAVQfK0iemoLkjOmAFImveiHYvBn6+4+kUkwW+K5atZHu7hFGR7uAXkZHu+juHmHVqo1ZFruQVDMl0jxRU8/iUMdnBcmZU7OkSHxDQ+MvBNULxdDQxOsvWHA1AwPtmHlnerNBBgbaWbDg6mwKXGCqmRJprqipZ1VhR81Sx+eQ4yQ75y4CvgJMA+4ys/9a8/pM4JvAe4CDwCfMbJ9z7kRgC7AY2Ghm19fbVtHHSda4xyLNoe9eevS/FcleNWhdvdpLKUujU3IW22i2ROMkO+emAbcBFwPzgU865+bXrPZp4CUzextwK3CLv/xXQA/wRzHLXghqlhRpPo05nh79b0WyFTX1LK64tdVFUbcm2Tm3BKiY2Qf85zcBmNmXAuts9dd51Dk3Hfg3oN38P+6cuxpYpJpk1biIiIhIMuvXe2kPwaB1cNBLPVuzpnHbUU1yfScDzwaeP+cvm3AdMzsEvAKcGKGA1zjnhp1zwyMjI2HfJiIiIlI6a9YcHax2daUTIG/eDMccU0mttjrPctFxz8zuNLNFZraovb292cVJlZolRUREJO+CHaX7+vrqdpQuojBB8n7glMDzOf6yCdfx0y1m43XgkxrKQxaRItDwUCLFlkVtdd6FCZKHgNOdc/OcczOAK4B7a9a5F1jlP74cGDAl3opIiyvTTW0w6K1+7qmCXg0PJdI64tzUatCB8EPAfRD4c7wh4L5hZl90zvUDw2Z2r3NuFvAt4BzgReAKM9vrv3cfcDwwA3gZeL+ZPTXZtorYcS+rBHsRaawydbQN5h92dzsGBqzuBAVl6NQjUgTB73dX19HP6ynyuTBpxz3M7H+Z2dvN7D+Z2Rf9ZWvN7F7/8a/M7ONm9jYzO7caIPuvzTWzE8zsN8xszlQBcqsJezelGhcRybvgDF7QF+oCWvbhoURaRRYz9BVRLjrutaqws+fp4BRpHWVtYqxUKnR3O0ZH+4G1jI7209099eceHPRqkHt6vN9l6vUu0mqS3NSWddCBUOkWWcp7ukUwdaLa/BA2dWLtWu/g7OnxBv8WkXwrchPjRKqtXKOj/bS1rQ2VahG3+VZEsqX0qIklTreQI3bu3Eh39wjOeUeWc110d4+wc+fGKd+nGhcRybNgkAu9dcdEDQ4PBZRyeCiRVpHVDH1FoyA5ok2brmZgoJ22tkGgj7a2QQYG2tm06epJ36ODU6Q1lamJMRj09vb21g16NTyUSOvQTW08SreIKUrqhEa3EBERkbKoVCot049jqnQLBckxVGuG58/fxlNPdSivR0RERMTXSv05lJPcQMHUiW3bOpQ6ISIiIlJACpIjUl6PiEwmj82LeSyTiBRPEYfPVLqFiEiD5LGJMY9lEpFia6XzjtItRERCWr/+6PSpwUFveSPWFxEpoiKeCxUki4gERJ1GPu7Y6WlavvwBnOuqafbsYvnyB5pWJhEptuC5s7e3t+65sxUoSBYRCYg6jXycsdPTdvPNy/1yeM2dAwNGW9sgN9+8vGllEpFiC547Dx+uFGIGTgXJIiI1urq8qVvXrfN+1zvJV9eHtaHWT1vwYgV9hbhYiUj+RT135p2CZBGRGlGnka+uv3TpttxMO5+3wF1Eii/quTPvFCSLiAREnUY+i7HT43SIKdrFSkTyLeq5sxUoSBYRCYg6FnoWY6dH7UxYxIuViORbEeeR0DjJIiItoBr4rl7t1QxPlWO8fr0XQAdfHxz0LlZr1mRTXhGRVpB4nGTn3EXOud3OuaedczdO8PpM59w9/uuPOefmBl67yV++2zn3gdifIm3bt8Ob3wz/9E/NLomIlN0E56MoHWLWrDn69a4uBcgikkM5jr/qBsnOuWnAbcDFwHzgk865+TWrfRp4yczeBtwK3OK/dz5wBbAAuAj4mv/38ueqq+CVV+C3f7vZJRGRkph0utYJzkd57BwoIpJYjuOvuukWzrklQMXMPuA/vwnAzL4UWGerv86jzrnpwL8B7cCNwXWD6022vczTLfzB9ieUs1QUESmWo6ZuneR8NEgnK9oG2bwZursdAwOmYd1EpLXlJP5Kmm5xMvBs4Plz/rIJ1zGzQ8ArwIkh34tz7hrn3LBzbnhkZCREkRropz+F3/zN8cvmzoUnnsi2HCIik5yPhj737cJ1iBGRkmuB+CsXo1uY2Z1mtsjMFrW3t2e78YUL4Y1vHL/sjW+Es87KthwiUgqVSsWfJjo4ZbTzUi8mOR+9Ovsv6O4e/57ubserr1YyLbuISMO0QPwVJkjeD5wSeD7HXzbhOn66xWzgYMj3Nt9LL8GCBXDPPd7vF19sdolEpKCOO67CwICNpVmYGQMDxnHHVbwVJjgfVSoVzMa/x8wmz2kWEWkFOY+/podYZwg43Tk3Dy/AvQKoza6+F1gFPApcDgyYmTnn7gW+7Zz778B/BE4HftyowjfM888feezN4yoikorqmMebN3vPg2MaAzofiUh55Px8VzdINrNDzrnrga3ANOAbZrbTOdcPDJvZvcDXgW85554GXsQLpPHX2ww8BRwCPmNmr6f0WUREcq+aT7xiBSxd+mDkDni9vb3pFlBERABNJiIi0hRr13pjHvf0eLPiiYhI9hJPJiIiIo1THfO4pweNeSwiklMKkkVEMhTMQe7vP5J6oUBZRCRfFCSLiGRoaAiNeSwi0gKUkywiIiIipaScZBERERGRCBQki4iIiIjUUJAsIiIiIlJDQbKIiIiISA0FySIiIiIiNXI3uoVzbgT4WZM23waMNmnbkj3t73LR/i4X7e/y0T4vl0bt7980s/aJXshdkNxMzrnhyYYBkeLR/i4X7e9y0f4uH+3zcslifyvdQkRERESkhoJkEREREZEaCpLHu7PZBZBMaX+Xi/Z3uWh/l4/2ebmkvr+VkywiIiIiUkM1ySIiIiIiNRQki4iIiIjUUJAMOOcucs7tds497Zy7sdnlkcZzzn3DOXfAObcjsOwE59wDzrl/8X+/pZlllMZxzp3inBt0zj3lnNvpnPusv1z7vICcc7Occz92zj3h7+8+f/k859xj/rn9HufcjGaXVRrHOTfNOfdT59zf+c+1vwvKObfPOfekc267c27YX5b6+bz0QbJzbhpwG3AxMB/4pHNufnNLJSnYCFxUs+xG4AdmdjrwA/+5FMMh4L+Y2XzgPOAz/vda+7yY/h3oNrOzgYXARc6584BbgFvN7G3AS8Cnm1dEScFngV2B59rfxdZlZgsDYyOnfj4vfZAMnAs8bWZ7zezXwHeAS5tcJmkwM3sIeLFm8aXAJv/xJuA/Z1kmSY+ZvWBmP/Ef/xzvQnoy2ueFZJ5f+E/f4P8Y0A1s8ZdrfxeIc24O8CHgLv+5Q/u7bFI/nytI9i6czwaeP+cvk+L7D2b2gv/434D/0MzCSDqcc3OBc4DH0D4vLL/pfTtwAHgA2AO8bGaH/FV0bi+WPwfWAIf95yei/V1kBtzvnHvcOXeNvyz18/n0Rv9BkVZkZuac03iIBeOc+w3gb4DPmdn/71U2ebTPi8XMXgcWOufeDHwPeGdzSyRpcc59GDhgZo875zqbXBzJxgVmtt8591bgAefc/xd8Ma3zuWqSYT9wSuD5HH+ZFN//ds6dBOD/PtDk8kgDOefegBcg/7WZ/b/+Yu3zgjOzl4FBYAnwZudctTJI5/biOB+4xDm3Dy9Fshv4CtrfhWVm+/3fB/Bugs8lg/O5gmQYAk73e8XOAK4A7m1ymSQb9wKr/MergP/RxLJIA/n5iV8HdpnZfw+8pH1eQM65dr8GGefcscByvDz0QeByfzXt74Iws5vMbI6ZzcW7Zg+Y2ZVofxeSc+6Nzrk3VR8D7wd2kMH5XDPuAc65D+LlN00DvmFmX2xuiaTRnHP/D9AJtAH/G+gFvg9sBk4FfgasMLPazn3SgpxzFwAPA09yJGfxZry8ZO3zgnHOnYXXcWcaXuXPZjPrd86dhlfTeALwU+AqM/v35pVUGs1Pt/gjM/uw9ncx+fv1e/7T6cC3zeyLzrkTSfl8riBZRERERKSG0i1ERERERGooSBYRERERqaEgWURERESkhoJkEREREZEaCpJFRERERGooSBYRERERqaEgWURERESkxv8FWmSRl1n8f04AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.subplots(figsize = (12, 4))\n",
    "mod.plot_response(sources, receivers, theoretical_times, marker='+', label='theoretical')\n",
    "mod.plot_response(sources, receivers, \n",
    "                  [[obs_times[i] for i in range(len(receivers)*j, len(receivers)*(j+1))] for j in range(len(sources))],\n",
    "                  marker='x', color='blue', label='picked times')\n",
    "# plot legend\n",
    "handles, labels = plt.gca().get_legend_handles_labels()\n",
    "by_label = OrderedDict(zip(labels, handles))\n",
    "plt.legend(by_label.values(), by_label.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "starting_times = model(starting_params, config)\n",
    "starting_times = [[starting_times[i] for i in range(len(receivers)*j, len(receivers)*(j+1))] for j in range(len(sources))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7faa31f5c640>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsEAAAD4CAYAAAANWzs4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAApYElEQVR4nO3de3xcdZ3/8ff3nMmlTdO0SWbSO+n9ltMWWgEF5U6roMVdROuquOKy8lMXvKO47Lrqiosr+FMX7AoIi1oBlbogFMQqW/khlosoVhc0E9J7ziTpJbe5fX9/zOlJQ0pvSTqTmdfz8cijmTPfnPOdnken7377me/HWGsFAAAAlBIn3xMAAAAATjRCMAAAAEoOIRgAAAAlhxAMAACAkkMIBgAAQMmJ5OOi9fX1trGxMR+XBgAAQAl5+umnfWtt9JXH8xKCGxsbtXnz5nxcGgAAACXEGNNyqOOUQwAAAKDkEIIBAABQcgjBAAAAKDmEYAAAAJQcQjAAAABKzrCEYGPMKmPMn4wxLxljrh2OcwIAAAAjZcgh2BjjSvqmpDdKWiRpjTFm0VDPCwAAAIyU4dgn+FRJL1lr/yJJxph1klZL+sMwnHvYXHPNNXruuefyPQ0AAICSs2zZMt188835nsYAwxGCp0pqPejxVkmnvXKQMeZKSVdK0owZM4bhssdue2ePOrtTebk2AABAqVqwOJPvKQxywjrGWWvXSlorSStWrLAn6roH3HzzzSq/77e6Z/PWE31pAACAknb9tefmewqDDMcH47ZJmn7Q42nBMQAAAKAgDUcI/o2kucaYmcaYcknvkPSTYTgvAAAAMCKGXA5hrU0bYz4kaYMkV9Lt1toXhjwzAAAAYIQMS02wtfankn46HOcCAAAARhod4wAAAFByCMEAAAAoOYRgAAAAlBxCMAAAAEoOIRgAAAAlhxAMAACAkkMIBgAAQMkhBAMAAKDkEIIBAABQcgjBAAAAKDmEYAAAAJQcQjAAAABKDiEYAAAAJYcQDAAAgJJDCAYAAEDJIQQDAACg5BCCAQAAUHIIwQAAACg5hGAAAACUHEIwAAAASg4hGAAAACWHEAwAAICSQwgGAABAySEEAwAAoOQQggEAAFByCMEAAAAoOYRgAAAAlBxCMAAAAEoOIRgAAAAlhxAMAACAkkMIBgAAQMkhBAMAAKDkEIIBAABQcgjBAAAAKDmEYAAAAJQcQjAAAABKDiEYAAAAJYcQDAAAgJIzpBBsjLnRGPNHY8zzxpgfG2MmDNO8AAAAgBEz1JXgRyU1WWuXSPpfSZ8e+pQAAACAkTWkEGytfcRamw4ePilp2tCnBAAAAIys4awJfp+kh17tSWPMlcaYzcaYzW1tbcN4WQAAAODYRI40wBjzM0mTDvHUddba9cGY6ySlJX331c5jrV0raa0krVixwh7XbAEAAIBhcMQQbK09/3DPG2PeK+liSedZawm3AAAAKHhHDMGHY4xZJemTks6y1nYPz5QAAACAkTXUmuBvSKqW9Kgx5jljzK3DMCcAAABgRA1pJdhaO2e4JgIAAACcKHSMAwAAQMkhBAMAAKDkEIIBAABQcgjBAAAAKDlD+mAcAAAA8EqZrk4l2+JKtbUo2RbXx7ffq+/d9Z18T2sAQjAAAACOSzbVq5T/slJtcSXbWpRui8smmtW7f284Jjouou7p5+ZxlodGCAYAAMBh2WxG6Y4d/au7flzGb1ZPxy4daBg8ttxRU9SRN8fIa6iQF3PlNTiKVTnSR27L8ysYjBAMAAAASZK1Vpn97eHKbsqPK+vHlUq0Kp1KSZIcI82tj8ibJHlLyuXFHHkNrmZNNHKMyfMrOHqEYAAAgBKU7etWym/Jhd22uNJ+XFm/WX3dXeGYyeMj8qJW3kxXXqxSXoOrhfWOxpSNnrD7agjBAAAARcxm0kq1b8ut7vq5wCu/WT2dbeGYcRWOlkYdeQtMEHYdeTFHdWOLdyMxQjAAAEARsNYqs8/vL2UIPqTW529TJpOWJLmO0fx6V95UyTu5Igi7rk6aMLpKGYYDIRgAAGCUyfbuD1Z1c1uQZdriyvjNSvb2hGOmT4jIi0pNsx15DZXyYq4W1DuqiJRW2H01hGAAAIACZdMppdpbw5XdVFuLlGhWz55EOKam0tWymJHX1F/K0BRzNaGSsHs4hGAAAIA8szar9J7dSgVhN9kWlxLN6k1sVzablSSVuUYLo668GUZerL+UYdp4I1NipQzDgRAMAABwAmV69g6o2834caX9FqX6esMxjRNzpQze/Ii8BldezNG8OkdlLmF3uBCCAQAARkA21adUojVc3U35cVm/Wb37OsMxtWNdeVEjb4kJ6nZzpQzVFYTdkUYIBgAAGAJrs0p37hywuqtEs3oSO8JuahURo0VRV96soJQh6KY2eRylDPlCCAYAADhKma6OIOjmdmXIJpqV8l9WOpmUJBkjza4NShkWlYWlDLNrHUUcwm4hIQQDAAC8QjbZO6ibmvWb1du1LxwTGxd0UzvZDbcgWxR1VFVO2B0NCMEAAKBk2WxG6Y7tYdgNd2Xo2B2WMowtd9QUdeTNNfIa+ksZYlXF202tFBCCAQBA0bPWKrO/vb9u148r68eV9F9WJp3rpuYYaW59RN4kyVtSLi/myGtwNWti6XVTKwWEYAAAUFSyfd2DShmyfrP6urvCMZPHB6UMM92gwYSrhfWOxpQRdksFIRgAAIxKNpNWqn1bbnXXzwVek2hWd0dbOGZchaOlUUfegv5ual7MUd1YShlKHSEYAAAUNGutMvvawl0ZUm1x2USzev2tymYykqSIYzS/3pU3VfJOrghLGWbUUMqAQyMEAwCAgpHt3a9kWzzcgizj576SvT3hmOkTgi3I5rjyYuXyGhzNr3NUESHs4ugRggEAwAln0yml2lsH7MpgEs3q2dMejqmpdLUsZuQ19ZcyNMVcTagk7GLoCMEAAGDEWJtVes/usIwh3IIssV3ZbFaSVOYaLYy68mYE3dQaHHkxV9PG000NI4cQDAAAhkWmZ++A1sEZP66036JUX284pnFiUMowPxJ2U5tX56jMJezixCIEAwCAY5JN9SmVaA1Xd1N+XPKb1bOvMxxTV+XKixp5Sx15sUo1xXKlDNUVhF0UBkIwAAA4JJvNKN25M/yQWiooZehp3xl2U6ssc7Qo6sibdaCUIbe6O2kcpQwobIRgAACgTFfHgA+pZf24Uv7LSqeSkiRjpNm1QSnD4rIw7M6pdeQ6hF2MPoRgAABKSDbZO6ibmvWb1du1LxwTGxd0UzvFlddQKS/malHUUVU5YRfFgxAMAEARstmM0h3bB6zuKtGs3o7dYSnD2HJHTVFH3lwjr6FCXsyV1+AoVkU3NRQ/QjAAAKOYtVaZ/e39uzIEK7t9fqsy6bQkyTHSvPoyeZOsmpaUh93UZk2kmxpKFyEYAIBRItvXHQTd/lKGrN+svu6ucMyU8UEpw8z+UoaFUUeVdFMDBiAEAwBQYGwmrVT71sG7MnT64ZjqCkfLYo68BSYIu46aYo7qxlLKABwNQjAAAHlirVVmX9vAul2/Wb2JbcpmMpKkiGM0v96VN03yTqkISxlOqmELMmAohiUEG2M+JukrkqLWWv9I4wEAKDXZ3v3Bqm5udTfj576SvT3hmOkTgi3I5rryYuXyGhwtqHdUTjc1YNgNOQQbY6ZLulDSy0OfDgAAo5tNp5Rqbx2wumsSzerZ0x6Oqal0tSxm5DUZebFKeQ25bmoTKgm7wIkyHCvBN0n6pKT1w3AuAABGBWuzSu/ZHbYOPrAFWU9iu2w2K0kqjxgtrHflnZTrptYUc+TFXE0bTykDkG9DCsHGmNWStllrf3ukP8zGmCslXSlJM2bMGMplAQA4oTLdewZ8SC3jx5VOtCjV1xeOmXmgm9qCSLjf7txaR2WUMgAF6Ygh2BjzM0mTDvHUdZI+o1wpxBFZa9dKWitJK1assMcwRwAATohsqk+pRGu4spvblSGu3n2d4Zi6Klde1Mhb6oSlDIujrqorCLvAaHLEEGytPf9Qx40xnqSZkg6sAk+T9Iwx5lRr7c5hnSUAAMPIZjNKd+4cvAVZ+86wm1plmaNFUUferFwpg9fgyos5mjSOUgagGBx3OYS19neSYgceG2PiklawOwQAoJBkujoGfEgt68eV8l9WOpWUJBkjzT5QyrC4LAy7c2oduQ5hFyhW7BMMACgK2WSvUn5LGHjTQfvg3q594ZiG6lzYbTrFCbupLYo6qion7AKlZthCsLW2cbjOBQDAq7HZjNLt23NlDH5L/xZkHbvDUoax5U6uqcS8gaUM0Sq6qQHIYSUYAFCQrLXK7E8MqNvN+nElE63KpNOSJMdI8+rL5E2y8pbkmkt4MVczJxo51O0COAxCMAAg77J93bkyBj8XdtNtcWUTcfV1d4VjpoyPyItaebPcsJRhYdRRZYSwC+DYEYIBACeMzaSVat86eFeGzv7PVFdXOFoWc+QtMEHYdeQ1uKodQ9gFMHwIwQCAYWetVWZf24BdGeQ3qzexTdlMRpIUcYzm17vypkneKRVh2D2phi3IAIw8QjAAYEgyvfuVaouHq7sZP/eV7O0Jx8yYEGxBNteV11AuL+Zofr2jcrqpAcgTQjAA4KjYdEqpRGtYxpDyW3KlDHvawzETxrg6OWbkNfWXMjTFXNVUEnYBFBZCMABgAGuzSu/ZfVDr4BbJb1ZP+3bZbFaSVB4xWljvyjsp2IIsKGWYWk0pA4DRgRAMACUs071nwIfUMn5c6USLUn194ZiZB7qpLYzIi7nyGhzNrXVURikDgFGMEAwAJSCb6lMq0XrQ6m5cSsTVu68zHFNX5cqLGnlLHXmxSnkNjhZHXVVXEHYBFB9CMAAUEZvNKN25c/AWZO07w25qlWWOFkUdebMGdlObNI5SBgClgxAMAKNUpqtjwBZk1o8r6b+sdCopSTJGmlNXJi9m5TWVhaUMsyc6ch3CLoDSRggGgAKXTfYo5b8cruym/Vzg7e3aF45pqA7qdpcfKGVwtSjqaGwZYRcADoUQDAAFwmYzSrdvD8Nu0o9Lfly9nbvDUoax5Y68qCNv3sBShmiVk9/JA8AoQwgGgBPMWqvM/sSAut2sH1cy0apMOi1Jcow0rz4ib7LkLS2X1+DIi7maOdHIoW4XAIaMEAwAIyjb152r2/WDUoZgG7JkT1c4Zsr4iLyolTfLDRpMuFoYdVQZIewCwEghBAPAMLCZtFLtWwfvytDph2OqKxydHHPUtLC/m5rX4Kp2DGEXAE40QjAAHANrrTJ725T0+3dlkN+s3sQ2ZTMZSVLEMVoQdeVNk7xTKsJShhk1bEEGAIWCEAwAryLTu1+pAx9SC8oYMn6Lkr094ZgZE4JdGea68hrK5cUcza93VE43NQAoaIRgACXPplNKJVrDMoZUsCtDz972cMyEMa5Ojhl5Tf2lDE0xVzWVhF0AGI0IwQBKhrVZpffsPqh1cIvkN6unfbtsNitJKo8YLax35TUGW5AFdbtTqyllAIBiQggGUJQy3XsGfEgt48eVTrQo1dcXjplZG5QyLIyE3dTm1joqo5QBAIoeIRjAqJZN9SmVaA1Xd9NtcdlEXL37OsMx9VVB2F3qyGuoVFPM0eKoq+oKwi4AlCpCMIBRwWYzSnfuDMsYwi3I2neG3dQqyxwtjjryZh8oZcit7jZUUcoAABiIEAygoFhrle3q7G8d3NaibKJZKb9V6VRSkmSMNKeuTF7MymsqC8Pu7ImOXIewCwA4MkIwgLzJJnuU8l8OA2/az7UP7uvaF45pqA5KGZY78mKV8hpcLYo6GltG2AUAHD9CMIARZ7MZpdu396/uHtiCrGNXOKaq3JEXdeTNC0oZGlx5MUfRKid/EwcAFC1CMIBhY61VZn9iwK4MWT+uZKJVmXRakuQ6RnPrXHmTJW9pfze1mRONHOp2AQAnCCEYwHHJ9nUr2dailB+UMgSBt6+nKxwztSYiL2rlzXbDUoYF9Y4qI4RdAEB+EYIBHJbNpJRq33ZQg4m4lIirp9MPx1RXODo55shbZIKwm+umVjuGsAsAKEyEYACSglKGvW0DdmWwiWb1JbYpm8lIkiKO0YKoK2+a5J3SX8owo4YtyAAAowshGChBmd79SoVhN9dNLeO3KNnbE46ZMSHYlWGuK6+hXF7M0fx6R+V0UwMAFAFCMFDEbDqlVKI1XN1NHdiVYW97OGbCGFcnx4y8JiOvoVJeLFfKUFNJ2AUAFC9CMFAErM0qvWf3QXW7LZLfrJ727bLZrCSpPGK0KOqqqfFANzVHXoOrqdWUMgAASg8hGBhlMt17+ld221qU8eNK+y1KJfvCMbNqy3K7MiyMhPvtzq1zFKGbGgAAkgjBQMHKpnqV8lvD1d20H5f14+rdvyccU1/lyosaecucsJRhcczVuHLCLgAAh0MIBvLMZjNKd+4csLqrRLN62nfKWitJqixztDjqyJtzoJTBldfgqKGKUgYAAI4HIRg4Qay1ynZ1DtiCLJtoVspvVTqVlCQZI82pK5MXs/KaysKwO3uiI5dSBgAAhs2QQ7Ax5sOSPigpI+lBa+0nhzwrYJTLJnuU8l8OA2/aD7qpde0LxzRUB1uQLXfCbmqLoo7GlhF2AQAYaUMKwcaYcyStlrTUWttnjIkNz7SA0cFmM0q3b+9f3T2wBVnHrnBMVbmjJTFH3jwjryFXytAUcxStcvI3cQAAStxQV4KvknSDtbZPkqy1u4c+JaDwWGuV2ZcIg26qrUU2EVef36pMOi1Jch2jefUReVOsvGX9W5A1TjByqNsFAKCgDDUEz5P0emPMFyX1Svq4tfY3Q58WkD/Zvi4l21rCut2M36ys36K+nq5wzNSaSG4LslluWMqwoN5RZYSwCwDAaHDEEGyM+ZmkSYd46rrg52slnS7pNZLuMcbMsgc+0j7wPFdKulKSZsyYMZQ5A8PCZlJKtW87qMFEXErE1dPph2PGV7paEjXyFpkg7Oa6qdWOIewCADCaHTEEW2vPf7XnjDFXSfpREHqfMsZkJdVLajvEedZKWitJK1asGBSSgZFirVVmb9uAXRlsoll9/lZlg25qEcdoQdSVN03yTqmQ1+DIi7maUcMWZAAAFKOhlkPcL+kcSRuNMfMklUvyD/sTwAjK9O4P9trNre5m/LgyfouSvT3hmBkTgl0Z5vZ3U5tf76jcJewCAFAqhhqCb5d0uzHm95KSki4/VCkEMNxsOqlUojVY3W1R6sCuDHvbwzETx7o6JWrkef2lDIujrmoqCbsAAJS6IYVga21S0ruGaS7AINZmle7cdVDdbq6bWm/7jrCUoTxitCjqymsMuqkFq7tTqillAAAAh0bHOBSMTPeeAa2DM36z0v7LSiX7wjGzaoNShoX9pQxz6xxF6KYGAACOASEYJ1w21auU3xqu7qb9uKwfV+/+PeGY+ipXXtTIW+bIa6iUF3O0OOZqXDlhFwAADB0hGCPGZjNKd+4csLqrRLN62nfqQOl4ZZmjxVFH3pyglCHmymtw1FBFKQMAABg5hGAMmbVW2a7OAVuQZRPNSvmtSqeSkiTHSHPqyuTFrLym8mC/XUezJzpyKWUAAAAnGCEYxySb7FGqrSUXeP0Wpf24sn5cfV37wjGTqoO63eUHShlcLYo6GlNG2AUAAIWBEIxDstnMQd3UWoJuas3q6dgdjqkqd7Qk5sibZ+Q19Jcy1I918jdxAACAo0AILnHWWmX2JXJh18/V7Vo/100tk0lLklzHaF59RN4UK29ZhbyYI6/BVeMEI4e6XQAAMAoRgktItq8rXNVNBluQZfy4kj3d4ZipNRF5UStvths0mHC1oN5RZYSwCwAAigchuAjZTOqgUoZ4UMoQV09nf0fr8ZWulkSNvEX93dS8mKuJYwi7AACg+BGCRzFrrTJ7dw9Y3VWiWb3+1rCbWplrtCAakTfdylueK2VoirmaUcMWZAAAoHQRgkeJTM++Aa2Dc93UWpTq6w3HnDQx2JVhXiT8kNq8OkflLmEXAADgYITgAmPTSaUSrWHYTbU1yybi6t3bEY6ZONbVKVEjb0l/KUNTzNX4CsIuAADA0SAE54m1WaU7dw1Y3VWiWb3tO8JShoqI0cKoK68x6KbW4MqLOZpSTSkDAADAUBCCT4BM954BrYNzpQwvK5XsC8fMqg1KGRZGwrA7t85RhG5qAAAAw44QPIyyqV6l/NZwdTftx2X9uHr37wnHRMcFYfdkR02xSnkxR4tjrsaVE3YBAABOFELwcbDZjNIdO8IyhqQfl/xm9XbskrVWkjSmzNHimCNvzsBShoZxdFMDAADIN0LwYVhrlenqCMoYcluQWb9ZyUSr0qmUJMkx0py6MnkNVp5XHuy362jWREcupQwAAAAFiRAcyCZ7cqu6bXGl/Bal2+LK+s3q694fjplUHZQyNDryGirlxVwtijoaU0bYBQAAGE1KLgTbbOagbmotQTe1ZvV07A7HVJU7WhJz5M038hoqwj1368dSygAAAFAMSiYEP/jgg7r7Ex/R7tZmZTJpSZLrGM2vj6hpipW3LNdNzWtw1TjByGELMgAAgKJVMiF43LhxWlbbJ29KfynDgnpHFRHCLgAAQKkpmRB81lln6azrL5aevTvfUwEAAECeUeQKAACAkkMIBgAAQMkhBAMAAKDkFExNcCqV0tatW9Xb2ztyF5n2Tin2lpE7f9GyqtzzF0175ssqS3bmezIAAABDVjAheOvWraqurlZjY6PMSG1P1tkidbePzLmLmLVWia5abdWnNPPJT+d7OgAAAENWMOUQvb29qqurG7kAjONmjFFdVUS9NbPyPRUAAIBhUTAhWBIBuIDl7g33BwAAFIeCCsEAAADAiVAwNcGv1Hjtg8N6vvgNFx32+XMuvVLXfui9Wnn268JjN//nd/WnP7folhs+ozvv+W994WvfliR99ur36/LL3nzU177xljv13R89JElKZzLa8mKz2p5/TLUTa9R42kWqHlcl13EUibja/NB3JUntHXv09quuVbx1uxqnT9E9t35ZEyeMP9aXHXry6ed127r1+s8b//G4zwEAAFAsWAkOrLlkpdat3zDg2Lr1G7TmkpVq79ijz920Vr9+4C499eB/6XM3rVVH596jPvcnrrpczz26Ts89uk5fuvZDOuv0U1Q7sSZ8fuO939Jzj64LA7Ak3fDNO3TemafqxV+t13lnnqobvnnHkF7fQxt/pVUHBXwAAIBSRggOXHrR+XrwsU1KJlOSpHjrdm3f5ev1p52iDb/8f7rg9aepdmKNJk4Yrwtef5oe/sUTx3Wd76/foDWXrDriuPUbfqnL33axJOnyt12s+x/+xaAx3/nBT3TJ+z6qC95xlRpPu0jfuGOdvvqtu3XyhWt0+sXvUXvHnnDsY5t+o/Nff6pe+NOfdepF79ayC96hJedfphf/8vJxvQ4AAIDRjBAcqJ1Yo1OXLdZDG38lKbcKfNmbL5AxRtt27tb0KZPCsdMmN2jbzt2DznH9jbfoJ4/88lWv0d3To4d/8YT++k3nhceMMbpwzQe1fNU7tfbuH4bHd/kJTW6ISpImxeq1y08c8py//9NL+tG3v6Lf/PRuXffl/9DYMZV69pHv67XLl+iu+x6QJPntHSqLRFQzvlq3/td9uvqKNbmV559+V9Mmx47hdwkAAKA4FGxNcD6suWSV1q3foNUrz9a69Rt0279ff0w//y+fuOqwz//3I4/rjBVLB5RCbPrx7Zo6OabdfrsueMdVWjCnUW84ffmAnzPGvOrOGee87jWqHlel6nFVqqkepzdf8AZJkrdwjp7/w4uSpEd++aQuPOt0SdJrly/RF//vbdq6Y7f+6o3nau6sGcf0GgEAAIoBK8EHWb3ybD226Sk987st6u7p1fIliyRJUyfF1Lp9Zzhu645dmjrp2FdQ1/3kkUGlEFODldhYfa3e+sZz9NRzL0iSGurrtGNXmyRpx642xepqD3nOivKy8HvHMaqoyD12jKN0JiNJeujnv9Kqc3L1wO986xv1kztu0pjKCr3p3R/Wzzc9dcyvAwAAYLQjBB9kXNVYnfO6FXrfRz+nNZesDI+vPOu1euTxJ9XRuVcdnXv1yONPauVZrz2mc+/Zu0+/fPJprV55dnisq7tH+/Z3hd8/8ssn1TR/tiTpLRe+QXfemytnuPPeB7R65VnH9ZqstXp+y4tatni+JOkvLVs166Rp+ocr1mj1yrP1/JYXj+u8AAAAo1nBlkMcaUuzkbLmklV66xUf07pbvhQeq51Yo3+85v16zUXvkiRd/5G/G1DScMD1N96iFUsX6S0XDg6sP35ooy58w+mqGjsmPLarLaG3XvExSbmt0955ySqtOucMSdK1H/xbXfaBT+m279+vk6ZN1j23fvm4Xs/Tz2/RyU0LwnKKe/77Uf3XDx9UWSSiSbE6febD7zuu8wIAAIxmxlp7/D9szDJJt0qqlJSW9H+stUf8//UVK1bYzZs3Dzi2ZcsWLVy48LjnclQ6W6Tu9pG9RoH5ws3f1pyZ0/WO1SuPPPgItrTs1sINlw3DrAAAQEn5yAtSzbS8XNoY87S1dsUrjw91JfjfJH3OWvuQMeZNweOzh3hODKPPXvP+fE8BAACg4Ay1JthKOtDGrEbS9iGeDwAAABhxQ10JvkbSBmPMV5QL1K/akswYc6WkKyVpxgy25QIAAED+HDEEG2N+JmnSIZ66TtJ5kj5irf2hMeYySbdJOv9Q57HWrpW0VsrVBB/3jAEAAIAhOmIIttYeMtRKkjHmLklXBw/vlfTtYZoXAAAAMGKGWhO8XdKB/cDOlcSmswAAACh4Q60J/jtJXzPGRCT1Kqj5HRb/PHgf3qGdb89hnz7n0it17Yfeq5Vn95c13/yf39Wf/tyiv3/XX+uqT/+r9u7vkus6uu7DV+jtx7DlWKK9U5de+Un95rcv6L2XvVnf+OK14XPfv/9h/evXb5cx0pSGqO7++hdUXztR7R179ParrlW8dbsap0/RPbd+WRMnjJe1Vldff6N++vNNGjumUt+56XM6xTv+reVSqZROu/hyPbPhe8d9DgAAgNFmSCvB1tpN1trl1tql1trTrLVPD9fETrQ1l6zUuvUbBhxbt36D1lyyUmPHVOqur31eL2y8Tw/f/U1d88//rs49+4763JWVFfr8J6/SV/7xIwOOp9NpXX39jdp477f0/M/u0ZKFc/WNO34gSbrhm3fovDNP1Yu/Wq/zzjxVN3zzDkm5FsgvNr+sFzet19ovf1ZXffpLg653LDY99ZzOeM3SIZ0DAABgtKFtcuDSi87Xg49tUjKZkiTFW7dr+y5frz/tFM2bfZLmzsrtaDFlUlSxuolqS3Qc9bmrxo7RmaeerMqK8gHHrbWy1qqru0fWWu3d16UpDVFJ0voNv9Tlb7tYknT52y7W/Q//Ijj+C73n0otljNHpy5eoc88+7djVNuC88dbtWvCGv9J7r/knzTvzEv3Nh67Tzx7/tc5Y/beae8ZqPfXs78OxD//iCb3xnDPU1d2ji979D1p6/tvVdO7b9INX/IMAAACgmBCCA7UTa3TqssV6aOOvJOVWgS978wVhu+EDnnr290qmUprdOLjrya133adb77rvqK9ZVlamW770GXnnvV1TTlmpP7z4F12x5hJJ0i4/oclBIJ4Uq9cuPyFJ2rZzt6ZPaQjPMW1yTNt2tg0690vxVn3s79+lPz7+I/3xpWZ97/6HtOn+2/WV6z+if/367eG4jU9s1tmvW66HNz6hKZOi+u3PfqDf//xerTrnVXe7AwAAGPUIwQdZc8mqsCTiQCnEwXbsatO7/+EfdcdX/1mOM/i37gPvuVQfeM+lR329VCqlW+66V89u+J62P7NBSxbO1Ze+fsegccaYQWH8SGZOnyJv4Vw5jqPF82brvDNPlTFG3oI5irfmepps27FbtRPGa+yYMfIWzNGjjz+pT33xa/qfXz+jmvHVx3Q9AACA0YQQfJDVK8/WY5ue0jO/26Lunl4tX7IofG7vvv266D1X64uf+qBOX75kWK733Av/K0ma3Thdxhhd9uYL9MTTv5UkNdTXhWUOO3a1KVZXK0maOimm1u27wnNs3bFbUydFB5274qDSC8dxVFFeHn6fzmQk5UohVp71WknSvNkn6ZmHvydvwRx99t/+Q/9y09pheY0AAACFiBB8kHFVY3XO61bofR/93IBV4GQypbde8TG959KLdOnFr7pt8jGbOimmP7zYHNYXP/r4r7VwzkxJ0lsufIPuvPcBSdKd9z6g1SvPCo6fpbvue0DWWj359POqGT8uLJs4Vg9vfEJvPPcMSdL2nW0aO6ZS7/rri/SJD7xHz/zuj0N9eQAAAAVrqFukjZwjbGl2XIwrOYd/yWveepHe+r5rtO7WG8Ox9zzwsB7/9bNKdO7Vd4Jg+p2bP69lTQsG/Oytd94jSfrA5ZcNOm/ja1Zp7/79SiZTun/DL/TI97+lRfNn658++gG94a/er7KyiE6aNlnfufkLkhPRtR/+O1329x/XbevW66Rpk3XPt74iORG96YKz9dONT2jOGas1dkyl7rjp84Nfk+MGvwbHjZEcJ/c4eC5jjV5q2aoF8+ZKkn73p7/oE5//qhzHUVkkoltu+Ozg8xpXGlt3+N9jAACAVzKFt+5qrD3xHYxXrFhhN2/ePODYli1btHDh8e93i2OzadMm3X333br11luP+me4RwAAYLQxxjxtrV3xyuOFuxKMEXXmmWfqzDPPzPc0AAAA8qLw1qYBAACAEVZQITgfpRk4OtwbAABQTAomBFdWViqRSBC2CpC1VolEQpWVlfmeCgAAwLAomJrgadOmaevWrWprG9z9DPlXWVmpadMGd8kDAAAYjQomBJeVlWnmzJn5ngYAAABKQMGUQwAAAAAnCiEYAAAAJYcQDAAAgJKTl45xxpg2SS0n/MLHrl6Sn+9JYMRwf4sX97Z4cW+LG/e3eOXz3p5krY2+8mBeQvBoYYzZfKg2eygO3N/ixb0tXtzb4sb9LV6FeG8phwAAAEDJIQQDAACg5BCCD29tvieAEcX9LV7c2+LFvS1u3N/iVXD3lppgAAAAlBxWggEAAFByCMEAAAAoOYTgV2GMWWWM+ZMx5iVjzLX5ng+OnzHmdmPMbmPM7w86VmuMedQY82Lw68R8zhHHxxgz3Riz0RjzB2PMC8aYq4Pj3N8iYIypNMY8ZYz5bXB/Pxccn2mM+XXw/vwDY0x5vueK42OMcY0xzxpjHggec2+LhDEmboz5nTHmOWPM5uBYQb03E4IPwRjjSvqmpDdKWiRpjTFmUX5nhSH4jqRVrzh2raTHrLVzJT0WPMbok5b0MWvtIkmnS/pg8GeV+1sc+iSda61dKmmZpFXGmNMlfVnSTdbaOZI6JF2RvyliiK6WtOWgx9zb4nKOtXbZQfsDF9R7MyH40E6V9JK19i/W2qSkdZJW53lOOE7W2scltb/i8GpJdwbf3ynpkhM5JwwPa+0Oa+0zwff7lPvLdKq4v0XB5uwPHpYFX1bSuZLuC45zf0cpY8w0SRdJ+nbw2Ih7W+wK6r2ZEHxoUyW1HvR4a3AMxaPBWrsj+H6npIZ8TgZDZ4xplHSypF+L+1s0gv8uf07SbkmPSvqzpE5rbToYwvvz6HWzpE9KygaP68S9LSZW0iPGmKeNMVcGxwrqvTmSz4sDhcBaa40x7BU4ihljxkn6oaRrrLV7cwtKOdzf0c1am5G0zBgzQdKPJS3I74wwHIwxF0vaba192hhzdp6ng5FxprV2mzEmJulRY8wfD36yEN6bWQk+tG2Sph/0eFpwDMVjlzFmsiQFv+7O83xwnIwxZcoF4O9aa38UHOb+FhlrbaekjZJeK2mCMebAIg7vz6PTGZLeYoyJK1dyeK6kr4l7WzSstduCX3cr9w/YU1Vg782E4EP7jaS5wadUyyW9Q9JP8jwnDK+fSLo8+P5ySevzOBccp6CG8DZJW6y1Xz3oKe5vETDGRIMVYBljxki6QLm6742SLg2GcX9HIWvtp62106y1jcr9Hftza+3fiHtbFIwxVcaY6gPfS7pQ0u9VYO/NdIx7FcaYNylXr+RKut1a+8X8zgjHyxjzfUlnS6qXtEvSP0m6X9I9kmZIapF0mbX2lR+eQ4Ezxpwp6X8k/U79dYWfUa4umPs7yhljlij34RlXuUWbe6y1/2KMmaXc6mGtpGclvcta25e/mWIognKIj1trL+beFofgPv44eBiR9D1r7ReNMXUqoPdmQjAAAABKDuUQAAAAKDmEYAAAAJQcQjAAAABKDiEYAAAAJYcQDAAAgJJDCAYAAEDJIQQDAACg5Px/ptJ7uTTS/m0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "starting_model=layered_earth_model(v=starting_params[0:2], z=[0., starting_params[2]], theta=[0., starting_params[3]], td=td)\n",
    "plt.figure(figsize=(12,4))\n",
    "starting_model.plot_model(model_span)\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "`gtol` termination condition is satisfied.\n",
      "Function evaluations 13, initial cost 5.6282e-02, final cost 6.1405e-04, first-order optimality 3.58e-09.\n",
      "Theoretical model:\n",
      "V1: 254 m/s, V2: 1487 m/s, z(0): -3.7 m, dip: -3°, td: 0.003 , -0.002 , 0.001\n",
      "Fitted model:\n",
      "V1: 229 m/s, V2: 1601 m/s, z(0): -3.5 m, dip: -3°, td: 0.001 , -0.005 , -0.001\n"
     ]
    }
   ],
   "source": [
    "res = least_squares(fun, starting_params, bounds=([150., 400., -10., -0.1, -0.005, -0.005, -0.005], [1000., 2500., -1., 0.1, 0.005, 0.005, 0.005]),\n",
    "                    args=(config, obs_times), verbose=1, loss='soft_l1')\n",
    "[v1, v2, z0, theta, td[0], td[1], td[2]] = [mod.v[0], mod.v[1], mod.z[1], mod.theta[1], mod.td[0], mod.td[1], mod.td[2]]\n",
    "print('Theoretical model:\\nV1: %.0f m/s, V2: %.0f m/s, z(0): %.1f m, dip: %.0f°, td: %.3f , %.3f , %.3f' %(v1, v2, z0, theta*180/np.pi, mod.td[0], mod.td[1], mod.td[2]))\n",
    "[v1, v2, z0, theta, td[0], td[1], td[2]] = res.x\n",
    "print('Fitted model:\\nV1: %.0f m/s, V2: %.0f m/s, z(0): %.1f m, dip: %.0f°, td: %.3f , %.3f , %.3f' %(v1, v2, z0, theta*180/np.pi, td[0], td[1], td[2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "calc_times = model(res.x, config)\n",
    "calc_times = [[calc_times[i] for i in range(len(receivers)*j, len(receivers)*(j+1))] for j in range(len(sources))]\n",
    "obs_times = [[obs_times[i] for i in range(len(receivers)*j, len(receivers)*(j+1))] for j in range(len(sources))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-2-13a568e85554>:83: UserWarning: marker is redundantly defined by the 'marker' keyword argument and the fmt string \".k\" (-> marker='.'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string \".k\" (-> color='k'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: marker is redundantly defined by the 'marker' keyword argument and the fmt string \".k\" (-> marker='.'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string \".k\" (-> color='k'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: marker is redundantly defined by the 'marker' keyword argument and the fmt string \".k\" (-> marker='.'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string \".k\" (-> color='k'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: marker is redundantly defined by the 'marker' keyword argument and the fmt string \".k\" (-> marker='.'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string \".k\" (-> color='k'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: marker is redundantly defined by the 'marker' keyword argument and the fmt string \".k\" (-> marker='.'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string \".k\" (-> color='k'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: marker is redundantly defined by the 'marker' keyword argument and the fmt string \".k\" (-> marker='.'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string \".k\" (-> color='k'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: marker is redundantly defined by the 'marker' keyword argument and the fmt string \".k\" (-> marker='.'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string \".k\" (-> color='k'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: marker is redundantly defined by the 'marker' keyword argument and the fmt string \".k\" (-> marker='.'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string \".k\" (-> color='k'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: marker is redundantly defined by the 'marker' keyword argument and the fmt string \".k\" (-> marker='.'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n",
      "<ipython-input-2-13a568e85554>:83: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string \".k\" (-> color='k'). The keyword argument will take precedence.\n",
      "  plt.plot(receivers, times[i], '.k', *args, **kwargs)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7faa31edcd30>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAD4CAYAAAAejHvMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAABS4ElEQVR4nO3df3xU5Znw/8+dhCSgiIq6j0oo/mqViKiQ1qlFZ8xKZaUCq6LFVgQpli26vro2yFODIhSa9MdjBR8X14LYXxZtUeq3rLVholhma4CVreK6q0AJsk9V1BakYUxyff84M8kkzK9zMufMmZnr/XrllczMmTlncmbOfZ3r3Pd1GxFBKaWUUkop1ass3xuglFJKKaWU32iQrJRSSimlVD8aJCullFJKKdWPBslKKaWUUkr1o0GyUkoppZRS/VTkewP6O+mkk2TUqFH53gyllFJKKVXktm3b9p6InJzsMd8FyaNGjWLr1q353gyllFJKKVXkjDF/TPWYdrdQSimllFKqHw2SlVJKKaWU6keDZKWUUkoppfrRIFkppZRSSql+NEhWSimllFKqHw2SlVJKKaVUH83NEA73vS8ctu4vFRok+5B+MJVSSimVT3V1MH16bzwSDlu36+ryu11e0iDZh/SDqZRSSql8CoVg3Tor/li0yPq9bp11f6nQINmH/PjB1Oy2UkopVVpCIZg3D5YssX6XUoAMGiT7lt8+mKWa3daTA6WUUqUqHIaHH4bGRut3//YwmUgkwvLly4lEIu5voMs0SPYpJx9MN/kxu+2FUj05UEop5Q2/JmPi7d26dXD//b0xQLp4JBKJUF9fT2NjI/X19QUfKGuQ7ENOPphe8Ft22wulenKglFLKG06SMV4E1m1tVntXXW1lhqurI6xbZ92fSmtrK9FolK6uLqLRKK2trbnboDzQINmH4h/MeCAWD9TSfTC94LfstldK8eRAKaWUN5wkY7y4ytnQYAXIiZnh6uoIDQ2pnxMMBqmsrKS8vJzKykqCwWDuNigPKvK9AepoyT6AoVB+g7PE7HZ8W0olq9r/5CDf+0IppVRxSUzGNDZmbmMSA+t586y2KZv2OBKJ0NraSjAYJBAIZNyuZJnhdM8LBAK0tLTYWoeviUjGH+Aq4A3gTeDuJI9XAT+PPf57YFTs/puAVxJ+uoEL061r3LhxovynqUlk06a+923aZN1fSOuwa9MmkZNO6t2u/reVUkqpgYq3LY2N9tqYxkYRsH5nsmXLFhk8eLCUl5fL4MGDZcuWLa48p9AAWyVFTJqxu4Uxphx4CJgEjAa+aIwZ3W+xW4EPRORs4P8ATbEA/CcicqGIXAh8GdgtIq84Ceb9pphGb2ajoeHoM9RQKHnW2yk/DpLza9cXu/w6MEQpVbz0uJMdp+OQwmF48MEIodByHnwwknb55mZYvbpvVnj16taM+yKeGV6yZAktLS2Fnxm2K1X0LL1Z4gDwXMLthcDCfss8BwRif1cA7wGm3zLLgG9lWl8hZJKdnllt2bJFli1bVpRnYsk4yQw7PZtW6WlGXCnlNT3uZMdpWzls2BapqrJikaqqwTJs2JaU/1u7y5cS0mSSswmSrwMeTbj9ZWBlv2VeBUYk3H4LOKnfMm8B56dYx1xgK7B15MiRnvxTBmLZsmVSXl4ugJSXl8uyZcsyPsfpZY5CDqqdHiDtXD5S2dMTEKWU1/S4446mJpE5c/rGInPmLMsqsA6FlmmAnCDvQTLwGeAPmdYlRZxJthtYF0tQbfcAqQdUd+kJiFLKa3rcyY7dNtxJnKD74mgDDZIH3N0Cq5/y/860LimQIFnE/Q+zF0G1k/fhRLZfSr005y49AVFKeU2PO9nxog3XfZHcQIPkCmAXcAZQCewAavst8zXgn2N/3wisS3isDHgbODPTuqQAguSBVGCw82G284VxctnF7jqcvAcRe19KJ/9bP1bE8CM9AVFKec2L446f2wA77aWTbpx2aBuQ2oCCZOv5/B3wX7FuFN+M3Xc/cE3s72rgSawScC8nBsRAEPi3bNYjBRAke/lBy/YL5rRDvtvZ6vj/ZsUK632sWLEl5/8r/eJnx88NiVKqOHlx3PFrG2C3vXS71FqplnHNxoCDZC9//B4ki/jzkoWTDvludwFparIC5MR1rFixJasMt1vZ6mJRqAcjpZTKNa/aALczw34cV2SHX09YMtEg2QV+7PxuZ5viQVbilzJTkOXXAYsiIrNmbRFYJrNmFebBxa5CPRgppZQb3G6T/ZYZ9qtCTFqlC5IzTiaijtZ/muJMBb/9uE3xiTs6OgIsXLiQjo5Axok7nBQVtzuPe7IpMDNZuTLCmjX1GNPImjX1rFyZeYKXQp8MJnFK0kWLSmeKcKWU6s9Jm2y3DbDTNjU3W21rYnvZ0REoiUlUEqfXnjevCNqkVNFzvn78nkn2YwbP6Tb58RKVk37PQ4Ysk7IyK1tdVlYuQ4YsS/lenHYBcdNAuk748YqGUkp5xUn7Z7edsdtu+DFO8EqxZZLzHhT3//F7kOzHvqDFFmTZCaqdHLzsBNVOtskuv5/kKKWUXznpOmi3G6CTQeileHwu1JMDDZJVUsX0JbYTxK5YsUVgsBhTLmAF1Zle22/1Kwv1YKSUUrnmRX9hJ+2lH5NQbvJjEjEbGiSro5R6kGVnoJ9X05Db2aZCPRgppVSueVVJwk7QW0xJqGKnQbI6SikHWXYPXk76r9md3MVudjtx2/xSMqiUP1NKqYFz0nVCxJtKEnbajVJPQhUaDZKVinF68LLbdcLO5C4D6Sft9oyJdmjDoJQaCLvHzkTZHtucnMzbPbZpwqCwaJCsVIxXBy87k7s4rbjhVQ1qN/tWK3doI60y8fIzYucYMmfOMgHruAbWVbhccnIyr9+n4qZBslJ5YHfQht0sr9szJjpZh0jpTeziR5rVd0+xBExz54oMG9b3MzJsmHV/LjkZVFdRMVigXCoq8t91QhU/DZLzrFgOqip7xVCDWsR+YG2nb7WT/of6Xcqe3c+g3f9tqe6LYjkB2bRJ5LjjrMC4sdH6fdxxue16JuKs3Fq2V+EGotQqT6jUNEjOs2I5qKrs+Hl/u5mtttu32kn/Qyf1SkuZk9H42X5u/fw5d1uxZCI3bRIZPNj6jAwenPtBzHaf49Vnqlj2n8oNDZJ9QL+UpaPYMmx2BsTY7VvtpP+hk0ogdieoKYaMqpNjjt3nlOJxLb6/E09A/LC/ndi0SaSqyuoeVVWVOWvrpMuWiLuD6uwq5ZM7lZwGyT6hl3dUKbDbBcRu/8Nly5bFAmTEmNwPWLSbrfZjozuQbbJ7nCq145rTbgpusxtgbtokcswxfU84jzkmfaDsRak1t/n1pFblT7oguQLliXAYHn4YGhut36GQ9aNUsQkEAgQCgayW7egIcMwxLVx8cSvbtwfp6Mj8vKFDg4hUYkwUkUqGDg2mXb61tZVoNEpXVxfRaJTW1ta02xcKwb33Rrj99vqedaxY0UIolPw5oRCsWwfTpkV63sf69YG8fr/b2qxtim9DfBvb2tIfd+wep0r1uGYMiFh/i1i3862uDqZPtz67Bw+2MnRokMWLA6xbl3z5J56Azs5WysqidHd3UVYWpbOzlSeeSP3ZDQQCtLS00NraSjAYzPp77icNDUffVyqfW+VAqug5Xz/FmEn2Y6ZJqXxz8r1w0ifZ6YBFu9lquxlxP00EI6J9krPl5+4WdgfOOik9qVSxQTPJ+eU0q6NUMXPyveh9TgCwsli1temf4yT7ZTdbvXp1K52dUaCLzs4oq1enz1ZHIhHq6+uJRqNUVlbS0tKS1XZFIhHXsnh294cfj2vNzVZGNXH94bC1TckyiE40NCTPoOfq9Qfi4MHW2Ge2C2OiHDzYSvx70p+1vQHGjSvszLBSrkoVPSf+AFcBbwBvAncnebwK+Hns8d8DoxIeuwCIAK8BfwCq062rGDPJShWiUu2756RPsp0qHU6mLRfx3wyLfuRFdtvLDLrd/ed0enulShkDySQbY8qBh4ArgX1AmzFmg4jsTFjsVuADETnbGHMj0ATcYIypAH4MfFlEdhhjhgMfDyysV0p5Id7HMZ4tDId7bxczu9nqtjZYvz5AdXVvRq6jI5By+bo6WLYsSEVFJRCloqKSJ58Msn59+u2y27faj9lqt8Wz2dOnw7x5VpY3MdudC15l0O3uv3AYFi8OsGJFS58+ybW1udkuL7L0SvlOquhZejPBAeC5hNsLgYX9lnkOCMT+rgDeAwzwd8CPM61DNJOslC8Veokvv2bD7U6Y4KT/qFczLPoxW+23ihtOJs6xu//c/qyXah90VfwYSAk44Drg0YTbXwZW9lvmVWBEwu23gJOAO4EfxYLo7UBDinXMBbYCW0eOHOnV/0UplQW/BRx2+LlhdzLRhxcDFu0EZk5LgrkVWDc1iXzve31P7L73vfzP4uhk4hw/llsr9JNmpZLJZ5B8F7A79vcQrL7J9enWp5lkpfyjGBpFP74HLyb6EHF3hkUR/2Wrv/c9EWOs38luJ+P0RMru/9bOxDlOMs9eKeSTZqf8ekXKbaXyvgcaJA+ku8WNwNqE5RqBb6RbnwbJSvmDn7OwdvmpYfdyog8n7E4G46dsdTyTnNiVJVMmWcT+CYjTQD/bMoF+/e758YTTC37dH24rlfc90CC5AtgFnAFUAjuA2n7LfA3459jfNwLrYn+fEOtmMST2Or8Frk63Pg2SlfKHYski+K1hd/p/9dv7iPNbttpJzWoRkauvXiUwUa6+elXGZZ1k0O32Q/fb/i6VgCkVv+0Pr5TC+x5QkGw9n78D/ivWjeKbsfvuB66J/V0NPIlVAu5l4MyE534Jq/zbq0BzpnVpkKyUypViadjtvg+/n+C4ma22060h7utfXxVb3vr5+tfTB8pOpzq3+zn00xUQv3+mvOCn/eGlbN93oX5GBhwke/mjQbJSKlcK9aDdn9334deTA6f7I9ugOp6tHTSoSowxMmhQVcas7aZNIoMGTewTJA8aNDHj/2rVqlUyceJEWbUqc+bZyfsuhQxeISnV/WHnffv1uJOJBslKKVVi/Nioz50rMmxY30Z02DDr/lyIl8qrrKwUY4xUVlZmLJXX1JQ8k5zuOW5X9SjUYKNYFcP+GMiJmp33bbdbkR+kC5LLUEopVXRCIWtCjSVLrN/5mio60Y03gghMmwaLFlm/Raz7c6GhwZqauaurCxGhq6uLgwdbM052MXnyXFatWsXEiRNZtWoVkyfPTbt8sold0mluhpUrrclBGhsbqa+vZ+XKCM3NyZePT1hSXR1h+fLlVFdHeiYsUd5LN4FMoYhPDhUOW7fjk0PV1aV+jpP3XV0d4aOP6gmHG/noo3qqqyO5exP5kCp6ztePZpKVUmrg/JhJFrG2Y/Bg6zrm4MG5L2XnRX9hJ+sYMmSZlJVZfaXLysplyJBlOV1H4vP8NrmL8gcvjglOxgTkG9rdQimlSoefLw9bAaPV+gwZ4l65NTuBopPgwe46VqzYIjBYjCkXsGZMTKWpyQo2EitozJmzLGO3EbuzMqrS42Twod0xAXYmzfEDDZKVUqqE+HXAYrwP8pAhVp/FIUO29OmjnIyTcmtOeFG5YNasLQLLZNas3AcbTrLVIpp5LiVOTwbt1CnftMmfk+Cko0GyUkqptLwIrOfOtQLkxOBvyJAtaQfueTE9sxeXoe2uw8kAKDvZahF/Tn2t3BH//NmZ2l7Eu5PUfEoXJOvAPaWUUo4G9gBEItbgskgk8wCds86CGTNa6ey0Br11dkaZMaOVs85K/ZxAIEBLSwtLliyhpaWFQCBg411lFn+f69bB/fdbvxP/D/laRygEd9wRIBxeyB13BLIaeDl/foBZs1oQWcKsWS3Mn5/+f2V3AKIqXG1tcO+9ERoarMGjDQ313HtvJOPgw2AwSGVlJeXl5VRWVhIMBj3ZXt9IFT3n60czyUoplR92s51O+wv7KXvpRQbdqzrJXuw/5S92uss4zQoXe5cctLuFUkqpbGTbb1ZEG123OK1P62Swpu6LwmX3JEdPipJLFyRX5DePrZRSyi9WroywZk09xkRZs6aSiy9Ofcm+uRmGDrUuxUajUSorKxk6NEhzMxnrEgcCgZx3mygm6erTpup24eQ5oPvCTyKRCK2trQSDwaz2SbLuMumeF++6ZGcdpc5YQbR/jB8/XrZu3ZrvzVBKqYJnp9ENh2Hy5OV0dDTS3d1FWVk51dVLePbZhUmDrHg/23vvjXDwYCtDhwZZvDjQJ1AbqOZmq0904uuFw1bglykQV6qQRCLWZDPxE85M/e+bm2HIEKuPcfw5zc0tHD4c0O+GTcaYbSIyPtljmklWSqkiZLfRbWuDpqYgDQ29meGmpmDKTGQ8Uzl9eoB58wJ8//vkNECG3sGE8ddNHACnlF/FT+6qq3tPUjs6AmlP7uxmha3vRoDm5pajTlJV7miQbJNmNpRShcBuo2sdvwKMG5f95djEqa8bG3M/9XVvIG6t5+GHcx+IK5VrdXUwbVqEjo56OjujVFRUUl3dwvr1qb9P8SoS8RPUTFUkvDhJVWgJOLuclklSSqmBslNuzWnppkAgwMKFC7PqrxgOW4FrY6P1O5dl0+ISA/F58zQIUP4XCsH117dy5Ih1knrkSJTrr29N+dltboaOjr6lDjs6AjQ3Z16PfjdclmpEX75+CqG6hReF55VSKpEX0zPb4dXU13q8VX5h5/u0ZcsWqagYLFAuFRXpv69Ov0v63cgNtARc7nkxhalSqi+/TrfshN0A1m8zX3mxL7wKxFX2SrVknN2TVLszJjqZkVG/G7mhQXKO6dmbUvlRLA1DMUzC4YViOikqBsX0GXTzJNXpccpO8k2/G7kz4CAZuAp4A3gTuDvJ41XAz2OP/x4YFbt/FPBX4JXYzz9nWpffg+RiaaSVKlR+PUnVma9UsfPb1Qyn3D5J9WqGRZUbAwqSgXLgLeBMoBLYAYzut8w/xANg4Ebg59IbJL+aaR1SQEGynr0plX9+6+6kM1+pUuDnz20hn6Rq8i2/0gXJ2VS3+DTwpojsEpEo8AQwpd8yU4C1sb+fAuqNMSaL1y44DQ1HjyANhbT8m1Je8aKiAtirJJGs3FoqTkayNzcf/T7DYTKOfi90pfq+/chpBQYvxGuCNzY2Ul9fn/Y7mzhTZLzyS3ymyEzsVH6xI91siSrPUkXP0pslvg54NOH2l4GV/ZZ5FRiRcPst4CSsTPJHwL8DLwATUqxjLrAV2Dpy5EhvTh2KiGa3VanwIuPS1CSyYkXfjNmKFVvSfp/sZNicvIdSzTSV6vv2Iy/3hRf9hVessNaxYsUW/UyVOAbY3WIgQXIVMDx23zigHTgu3fr83t3Cj7QhUaXC6QmhnUZ30yaRIUOWSVmZ1eiWlZXLkCHLUn6f4tuUuA43+h+Wap/FUn3ffuT2vnBygirirPKEfqZU3ECD5ADwXMLthcDCfss8BwRif1cA7wEmyWu1AuPTrU+DZGf0S69Uck76Ua5YsUVgsBhTLmA11Kl4MZJ9IM8pBqX6vv3ISQWGbE8g7Z6gJrKbfdbPlIobaJBcAewCzqB34F5tv2W+Rt+Be+tif58MlMf+PhN4Gzgx3fo0SHZOv/RKHc3pIJ1Zs7YILJNZs7LLPjupcaqZ5MxK9X37kZPP+bBhW6SqyjpJraoanLFmsJ0TVKf0M6USDShItp7P3wH/FetG8c3YffcD18T+rgaexCoB9zJwZuz+a4HXsMq/bQe+kGldGiQ7o196VUrsznxlN5Ps5PuU7Umq9knOXqm+bz9yui/mzFkmYJ2kQrnMmZP5JNXOCapd+plS/Q04SPbyR4Nk+/RLr0qJ29MzDySIzSaodtKvulQH55bq+/YjJ33vRexNzyziTb9n/UypRBokFzn90qtC5rfpme1+n/QkVZUKp1dlsp2e2cl3Sds/NVDpguRs6iQrn9PazcoNXtSotVPfNC4Y7FvjNBgMplzWyXuw+33SGqeqkLlVDxys79r06bB0KVx5pfV7+vTUtc2dfJfq6vq+ZnyddXUZ345SmaWKnvP1o5lkpXJvINOk2s2Q+mnmK83yKpWa3amW7ZZnmzRJ5Pbb+z7n9tu3yKRJuX0fOiZHDQTa3UIpS6lemnMaLNptfPw4PbM2oKoU2D22NTVZg+oST1LnzFmWsVuRnUk4BlLSzS6t7qSc0iBZqRgvMot+DcSdBot2Rpo7yQzb7ZPshDagqtjZDWKdlGdzcgzRkm7K7zRIdoEXDbtyh9sHVD9f4rcbLNpt4LzIDNulDagqFXa/r3YG1cU5OeHUkm7KzzRIzjGngUAhB9Z+zY465XZm0Y+Bmd1tsnup1GmJKDdpA6pKybJly2IBMmJMdldy7BwLnRzXtKSb8jsNknPM6SVlN2u7uq2Yuil4FcD66RK/k/6Edgfq+DEg1QZUlRInmeRsj4VOvt9+PCYo1V+6ILkiLyU1Cly8BFU0Gs1YgiouWemcQCCQcvl4aaz4OlpaWtIuH39Oa2srwWAw47J2xUvxTJ8O8+bBww/3LdWTC/FSPvHXjZfyWbcud+tIfM1QyPpJvJ3L9Tz8MDQ2Wr/j68qXtja4994IDQ29n6nm5hba2gIpt8sqeRZg3LiWrD5XXnxG7EpWti3f+0IpN4TDsHhxgBUrWjh4sJWhQ4MsXhygtjb5593usTBdebZU3ycnzykkH3/8Mfv27aOjoyPfm6KyUF1dzYgRIxg0aFD2T0oVPefrpxAyySL2s7x2M8l2s9VedQEp9G4KXmQWvcyeeFFuzS4/ZdCVKhVOqlvoVZaB2bVrl7z77rvS3d2d701RGXR3d8u7774ru3btOuoxtLuFP9gJaNwOqp2sw8kgDycKPcjyqvEp9HJr2kgrpQrZzp07NUAuIN3d3bJz586j7k8XJOuMex4KBAIsXLgwY1eI5mbo6AjQ0tLCkiVLaGlpoaMjkHaWMDuzkMXZmT0pHIZp0yJ0dNTz4ouNdHTUM21aJOXMSYnszOjUv5tCNq/vJi9mbIuz838C+7NfBQJ9P1O57pKTePn2/vt7u16k2oc6U5ZSqtAZY/K9CSpLjvZVqug5Xz/FnEnOlhczndkdlGW38LyTdTgdXOZmNtKrrhNOB3b6qdzaQGb181MVEKVUfvlp0Ho6ybKSyt/sZpLzHhT3/9Eg2eJVLV87AamT7hl2S4jZnfbUyfuwy+m+yPZA7/QExG/l1pwq9O41Sqnc8dvJfzp+DpJvvfVWee2111I+fu+998p3vvMdR6+9Zs0a+drXvnbU/eFwWH73u9/13H744Ydl7dq1jtbhFg2Si4gfB8nZPcO3W5LISd9qJzM62QlgN23quy+yCUbtHOgHMvNVoZdW0kyyUiqRVwOMc8FOkOy3MRhuBMkDeU2vaJBcJIqplq+d2ZacZBHsFtC3H8CKVFWtkjPPnChVVatk2LDMXUDsZoadDIos9ACzWAJ9pdxSKN0OcsXJ1cR8shMku3G82717t3zqU5+SGTNmyLnnnivXXnutfPTRRyIicvnll0tbW5uIiGzcuFEuuugiueCCC+SKK64Qkb4B7SOPPCJXXXWVHD58WH70ox9JXV2djB07VubOnSudnZ0iIrJ69Wo555xzpK6uTubMmXNUkLx79275m7/5GznttNNk7Nix8uKLL/ZZx+WXXy533nmnjBs3Ts4991x5+eWXZdq0aXL22WfLN7/5zZ7XSbb+zs5OmTlzptTW1sr5558v3//+9x3/zzRILgJeBQ9eBFmFnq3etMkKkIGen6qqVRmzvMOGbZFBgyrFGCODBlWmDXydZqtFCrurgt8yK0r5RaEFi7niRfe5XLLb3SLXbe7u3bsFkJdeeklERGbNmtUnKG1ra5N33nlHRowY0VP67MCBAyLSGySvWLFCrrnmGuno6JCdO3fK5MmTJRqNiojIvHnzZO3atbJ//36pqamRd955R44cOSKf/exns8ok9w+SGxoaRETkgQcekFNPPVX2798vHR0dcvrpp8t7772Xcv1bt26Vv/3bv+153Q8++MDx/8yV6hbGmKuMMW8YY940xtyd5PEqY8zPY4//3hgzqt/jI40xh4wxd9kfWlh60hVgzxW7lQi8XEe2VUAAbrsN7rnHKqD/rW8tYcWKFu65J8Btt6V+jp1KIG1tMGbML/rcN2bML9Lui1AIli6Fjz82iFi/ly5NXTy/rg6mTYPvfS9CKLSc730vwrRpmas8+K0SiF1Oq4AoVezq6mDBglaOHLGq1xw5EmXBgtair/wSb+sWLw7w178uZPHiQN4nJMqlUMiaaGnJEut3Lt5XTU0Nl156KQBf+tKXeOmll/o8/m//9m9cdtllnHHGGQCceOKJPY89/vjjbNy4kaeeeoqqqipaWlrYtm0bdXV1XHjhhbS0tLBr1y5+//vfEwwGOfnkk6msrOSGG25wtK3XXHMNAGPGjKG2tpZTTz2VqqoqzjzzTNrb21Ou/8wzz2TXrl3cfvvt/Ou//ivHHXeco/U7kTFINsaUAw8Bk4DRwBeNMaP7LXYr8IGInA38H6Cp3+PfBzYOfHNLgxfBgxeBuBfrABCB2lorsK6tDSCSfvlAIMADDzxAfX09DzzwQNpgvKEBLrvs2j73XXbZtRn3xcGDrRjTCQjGdHLwYGva5T/+OMLhw/W0tjZy+HA9H3+cvgycFyc5Sqn8CIWgqSlId3clxpTT3V1JU1OwaILFdNwIJP3CjcRG/7JmdsqcjRkzhj179rBv3z7A6lkwc+ZMXnnlFV555RXeeOMN7rvvvoFvZExVVRUAZWVlPX/Hb3d2dqZc/wknnMCOHTsIBoP88z//M3PmzMnZNmWSTSb508CbIrJLRKLAE8CUfstMAdbG/n4KqDexPWWMmQrsBl7LyRarnPAiEPdiHatWwdNPWwHiokXW76eftu5PJRKJcOedd9LS0sKdd96Zti5xOAyPPz6Xr399FRMnTuTrX1/F44/PzXhwGzo0iIjVwIlUMnRoMOWybW0wY0YrxkQR6cKYKDNmtKY9mfDqBEQplR/z5weYNasFkSXMmtXC/Pm5rWvuV4V+hSwVtxIbe/fu7WnDfvrTn/K5z32uz+OXXHIJL774Irt37wbg/fff73nsoosuYtWqVVxzzTXs37+f+vp6nnrqKd55552eZf/4xz/ymc98hhdeeIEDBw7w8ccf8+STTybdlqFDh3Lw4EHH7yXV+t977z26u7u59tprWbp0Kdu3b3e8DrsqsljmdKA94fY+4DOplhGRTmPMn4HhxpgOYAFwJZCyq4UxZi4wF2DkyJFZb7xSYAWIX/hChCVLWpk1K0golL4xSTYJR6pscm8w2vMRZfJk6/5UGY5w2LpcuGJFCwcPtjJ0aJDFiwPU1iZ/TkMDrFwZD6qjiFQydmyQ+fNTv4dkJxqhUHFlXZQqFs3NVheKxO9nOGwdR1IlDcJh+NWvAjQ2BnqCxWL/ficGkvHjWeLtQpYusTGQ9/apT32Khx56iNmzZzN69GjmzZvX5/GTTz6ZRx55hL//+7+nu7ubU045heeff77n8c997nN897vf5eqrr+b5559n6dKlTJw4ke7ubgYNGsRDDz3EJZdcwn333UcgEOD444/nwgsvTLotX/jCF7juuut45plnWLFihe33Mnr06KTrHzx4MLNmzaK7uxuA5cuX235tx1J1Vo7/ANcBjybc/jKwst8yrwIjEm6/BZwEfBeYHrvvPuCuTOvTgXvKLrsD99yuw2l3QFqhDVZRStljdzB2qVZ+KbTBvPmuk7x7926pra3N6zYUGrsD97LJJL8N1CTcHhG7L9ky+4wxFcAw4ABWxvk6Y0wzcDzQbYzpEJGVtiJ5pVIIh60BLmVlUbq7uygrswa41NYGUp6dx6dnbm1tJRgM5nx6ZrtZ3t4MQwCwtqW2duAZBqWUP8SzhtOnW/1sH344fXbUrayj3+kVMuU3RjKMcooFvf8F1GMFw23ADBF5LWGZrwFjROSrxpgbgb8Xken9Xuc+4JCIfDfd+saPHy9bt2518l481d7ezp49exg1ahQ1NTWZn6Bc0dwMQ4ZEaGioJxqNUllZSXNzC4cPB5IecJ1c9lRKqVxYtMgakNbYaPVLVYXt9ddf57zzzsv3ZvhONBrlyJEjVFVVUVlZme/N6SPZPjPGbBOR8cmWz5hJFquP8XzgOaAcWC0irxlj7sdKUW8Afgj8yBjzJvA+cOMA34evtbe38/jjj9PV1UV5eTk333yzBsp5YgW2AcaNyy4zXFfXt49bYh84pZTKxOmJdv8BaZohza1IJOLa1cFSZyfojUajHDhwABHh0KFDDB8+3HeBsh3ZdLdARH4N/LrffYsS/u4Ars/wGvc52D5f2rNnD11dXYgIXV1d7NmzJ6sgWbPP7gkEAlkdGO1e9lRKqUROTrSLeUCaH0QiEerre68mtrS0aKCcgt0sr92g98iRI/GxaYgIR44cKeggOavJRFRfo0aNory8HGMM5eXljBo1KuNz4tnncDjM448/Tnt7e8bnKHcUUx3OSCTC8uXL05axU0rlTuKJdrzsZKZgN97HuLra+r5WV0e0ZGMOJatYVCqi0SgHDx4kGo1mteyBAwc4ePAgBw4cyOo5yYLedKqqqnpqNRtj+tRDLkRZZZJVXzU1Ndx88822ssJOss+aeXZHsVz21OyJUvmReKLd2Jj5+NHQkPz72tCg39eBam626tJXVlb2/G+HDg3S3Fx440z8mOWtqqri0KFDiEhWQW9lZSXDhw/3bZ9kuzST7FBNTQ0TJkzIOni1m312knlub29n8+bNBZulbm62AtjE7Gg4bN2fK8UyU11zM6xe3Td7snp1a07/V0qp5JxMeFHK2U431dVZdembm1tYsmQJzc0tLF4c8MUU3vnK8j7wwAMcPnz4qOUzZXkfe+wx9u/f33N7zpw5vPnmmwwfPpyhQ4dm3b/YOlEZ6lqAPGrUKN57770BL5MNzSR7xG722W7muRgGE9bVwbRpETo66unsjFJRUUl1dQvr1+cu21IspZXq6mDZsiAVFZWA9b968skg69fne8uUKm5O+xcHg32zncFg0LNtLma93V8CzJsX4Pvf90df73xmeR944AG+9KUvMWTIkD7Lp8vydnV18dhjj3H++edz2mmnAfDoo4/2eW4p0kyyB+IZ0sTsc6YMqd3Mc7KgOht+yj6HQnD99a0cOWJlW44ciXL99a05Pdh5MVW2F0IhWL8+QHV1C5ddtqTnZCLfDYNSxc7plPDx+uxLlizRrlE55tU4EzvtZaosb6rXcNKX9+OPP2b27Nl8/vOf52//9m9Zv349Dz74IPv37ycUChGK/SPmzZvH+PHjqa2t5Vvf+lZPlnfUqFEsWLCAiy++mJ/97Gds3bqVm266iQsvvJC//vWvBINB4iV5jz32WL75zW8yduxYLrnkEv70pz8B8NZbb3HJJZcwZswY7rnnHo499tijtnPPnj2ce+653HLLLXzyk5/kpptu4re//S2XXnop55xzDi+//DJgTUM9depULrjgAi655BL+4z/+A4ADBw4wceJEamtrmTNnTs//FeDHP/4xn/70p7nwwgu57bbb6Orqyvh/syXVLCP5+inGGfeczp60d+9eefHFF2Xv3r0Z17F3715ZunSpLF68WJYuXerac9y2ZcsWqagYLFAuFRW5nw2v2DQ2ioD1WymlSlG8TW1szH5mQjvta3z5/u1luhn3jhw5Ivv375e3335b9u/fL0eOHMnY5h45ckT+8pe/yJEjR7LapqeeekrmzJnTc/vDDz8UEZFPfOIT8u677/bcf+DAARER6ezslMsvv1x27NjRs1xTwnSGl19+ubS1tSW9DciGDRtEROQb3/iGLFmyRERErr76avnpT38qIiIPP/ywHHPMMUdt5+7du6W8vFz+4z/+Q7q6uuTiiy+WWbNmSXd3tzz99NMyZcoUERGZP3++3HfffSIi0tLSImPHjhURkdtvv10WL14sIiLPPvusAPLuu+/Kzp07ZfLkyRKNRkVEZN68ebJ27dqk/4M4N2bcUwPktOxYTU1N1l0mimUwYUdHgGOOaeHii1vZvj1IR4dmW1IplgGISinlVLz7y6OPtnPiiXsYN24U06fXpGxjm5vhU59q59VXe7snnn/+zbzxRk3aK4rJ2suTTjop5fLJujZkanMrKyttdWsYM2YM//RP/8SCBQuYPHkyEyZMSLrcunXreOSRR+js7OR//ud/2LlzJxdccAEAN9xwQ1brqqysZPLkyQCMGzeO559/HrDGED399NMAzJgxg7vuuivp88844wzGjBkDQG1tLfX19RhjGDNmTM+V75deeolf/OIXAFxxxRUcOHCAv/zlL7z44ov88pe/BODqq6/mhBNOAKClpYVt27ZRF+uA/te//pVTTjklq/eTLQ2SPWJ3NLQTdoJq6O3SET9QZDuY0K1+z/GDndVtIHBU3z/VS+uuKqWKlZ1kTFubFSAnBr2PPnozbW01SY+FdXXQ1LSHQKALEDo7u3j44T0sWJB+Pcnay0OHDqV9Tv+g126bm8knP/lJtm/fzq9//Wvuuece6uvrWbRoUZ9ldu/ezXe/+13a2to44YQTuOWWW+jo6Oh5/JhjjslqXYMGDerpDlJeXk5nZ6etbU3sPlJWVtZzu6yszPZrxYkIM2fOZPny5Y6enw3tk+wRJ6Oh3RbPPodCoawCXrf7PTvt61eK9H+llCoUdvry2q3s1NAAJ57Yt2068cQ9KbPCVsJqFB9/XI6I4eOPy5k3b1RWV3bttJduvUai/fv3M2TIEL70pS/xjW98g+3btwMwdOhQDh48CMBf/vIXjjnmGIYNG8af/vQnNm7cmPL1Ep+XrUsuuaQn+/vEE084fCeWCRMm8JOf/ASwqsGcdNJJHHfccVx22WX89Kc/BWDjxo188MEHANTX1/PUU0/xzjvvAFaf5j/+8Y8D2ob+NJPsAT9n/exkn52cBdvJPscPaolZhFAoeTag1CVrALS7hVLKbXa73Nm9AumkG6DdtmnKlBqeeaaWTZs2csUVk5gyJftujQMNbHPxGnF/+MMf+MY3vkFZWRmDBg3i4YcfBmDu3LlcddVVnHbaaYTDYS666CLOPfdcampquPTSS1O+3i233MJXv/pVBg8enPUEVfFKGt/61re46qqrGDZsmOP3c9999zF79mwuuOAChgwZwtq1awG49957+eIXv0htbS2f/exnGTlyJACjR49m6dKlTJw4ke7ubgYNGsRDDz3EJz7xCcfb0J+RhFGCfjB+/HiJj6YsFs3N1iWexAAmHLayfoVWVcHuAXLz5s2Ew+GeEjWhUChlv6n46xd6KTullCpGTo7PXrUBdtqmlSsj3H57PcZEEalkxYoW5s+3P/7l9ddf57zzzrP9vGJy+PBhBg8ejDGGJ554gp/97Gc888wz+d6slJLtM2PMNhEZn2x5zSR7oJiyfm73e3aSRQCdnVAppZywc+z0IsvrZBB6/HnZLBsOw4IFrZSVRenu7qKsLMqCBa3U1moJTSe2bdvG/PnzERGOP/54Vq9ene9NyikNkpWr7B7w3O7SUeoikQitra0Eg0Gt06rUAPj1u2Qn6LV77HRyfHYS9OayS0J/bW3Q1BSkoaF3YpempmDBTSjlFxMmTGDHjh353gzXaJCsXLdv3z5eeuklKioqMh74iqWUnd80N8OQIREaGup7Gobm5hYOHw4UXJcfpfItEolQX9/7XXJrchC/9f91O8vrBet4F2DcuBZfnuQof9EgWbnKSWNS6KXs/KiuDiZPtmYz7O62ZjNcsKCVZ5/VxkEpu1pbW4lGrZlBo9Eora2tWQVabmZ5wX7Q6zQzXAzHy0AgoMGxykiDZOUqp42JHXazG6WYeQ6FrEuMt99eiTFRurutS4x6eVEp+4LBIIMGDUJEGDRoEMFgMONz/FjlwWlmWKlSoUGyyorTCh3BYJDKyt6+X9k0Jk64WcrOi9HWXpg/P8D27S2sWdPKrFlBR6O5lSpWdr6vI0aM4Oabb+att97irLPOYsSIERlf36ssr5/6/ypV6LKaTMQYc5Ux5g1jzJvGmLuTPF5ljPl57PHfG2NGxe7/tDHmldjPDmPMtBxvv/JIXZ1V2zk+CUq89nNsNsiUAoEALS0tLFmyxLV+e3Z5MYmK3YL4XgiH4Ve/CtDYuJBf/SrgiwltlPIDu9/XPXv2cPrppzNhwgROP/30rI4J8aDXGGMry2t34omamhomTJiggW+JePDBBznvvPO46aab2LBhA9/+9rcBePrpp9m5c2fPco899hj79++39dp79uzh/PPPz+n2JnPsscfmZBk3ZMwkG2PKgYeAK4F9QJsxZoOI7ExY7FbgAxE52xhzI9AE3AC8CowXkU5jzKnADmPMr0TE2RyEKm/iM7pNn25Nr/3ww5knQ+nNPvf2/fJLfWi3J1HxW5cOP09oo1Su2f0uaZZXFar/+3//L7/97W97rmZcc801gBUkT548mdGjRwNWkHz++edz2mmn5W1bC1E23S0+DbwpIrsAjDFPAFOAxCB5CnBf7O+ngJXGGCMihxOWqQb8NXOJssWazhOWLLGm184UXMWzz/FALDFQKyROGje/DSZMN421BsnK79we8Oa3Wr6qeOWydOBXv/pVdu3axaRJk5g9ezYnnHACW7duZcaMGWzYsIEXXniBpUuX8sUvfpGtW7dy00039cymt3PnTr7+9a9z6NAhTjrpJB577DFOPfVUtm3bxuzZswGYOHFi0vW2trZy7733cvzxx/OHP/yB6dOnM2bMGH7wgx/w17/+laeffpqzzjqLPXv2MHv2bN577z1OPvlk1qxZw8iRI9m9ezczZszg0KFDTJkypc9rf+c732HdunUcOXKEadOmsXjx4gH9jwZMRNL+ANcBjybc/jKwst8yrwIjEm6/BZwU+/szwGvAIWBainXMBbYCW0eOHCnKnzZtEjnpJJHGRuv3pk3uPKdY7N27V1588UXZu3dvxmVffPFFWbx4sdx3332yePFiefHFF3O+DqX8ws7ndu/evbJ06VJZvHixLF26NONz9Luk3LJlyxZZtmyZbNmyRUREdu7cafv5gwcPlvLychk8eHDP6wzEJz7xCXn33XdFRGTNmjXyta99TUREZs6cKU8++WTPcpdffrm0tbWJiEg0GpVAICDvvPOOiIg88cQTMmvWLBERGTNmjLzwwgsiInLXXXdJbW3tUesMh8MybNgw2b9/v3R0dMhpp50mixYtEhGRBx54QP7xH/9RREQmT54sjz32mIiI/PCHP5QpU6aIiMgXvvAFWbt2rYiIrFy5Uo455hgREXnuuefkK1/5inR3d0tXV5dcffXVPdsSX2agku0zYKukiIFdH7gnIr8Hao0x5wFrjTEbRaSj3zKPAI+ANS2129uk7HN6ud5u9rmYuN2lw0nGzG+DCVXh81stXyffJdAsr0ovWTnT448/3tZreFHtKRtvvPEGr776KldeeSUAXV1dnHrqqXz44Yd8+OGHXHbZZQB8+ctfZuPGjUlfo66ujlNPPRWAs846qyfrPGbMGMKxAS+RSIRf/vKXPa/VEOtn+bvf/Y5f/OIXPfcvWLAAgN/85jf85je/4aKLLgLg0KFD/Pd//3fP9uRDNkHy20DikWNE7L5ky+wzxlQAw4ADiQuIyOvGmEPA+VhZY1VAnF6uD4et/suNjdbvQp2O221eTKJSivWhlX2FXstXy5qpXGtuhv/+774B7urVrdx661Rbr+NVtadMRITa2loikUif+z/88MOsX6Oqqqrn77Kysp7bZWVldHZmHnZmjEm6XQsXLuS2227Lejvclk11izbgHGPMGcaYSuBGYEO/ZTYAM2N/XwdsEhGJPacCwBjzCeBcYE9Otlx5qqHh6OA2FEo/AC8x+3z//b0D/7SqQnJ2R6XbHS3vpEoHWANAbrvtNp5++umsllf+0d7ezubNm7OuruKkyoPdz5QXVR6cVHiIRCIsX778qMBBqbo6ePLJIBUVlZSXl1NRUcmTTwZJiBOz4mW1p6FDh3Lw4MGktz/1qU/x7rvv9nzWP/74Y1577TWOP/54jj/+eF566SUAfvKTnwxoGz772c/yxBNP9LzWhAkTALj00kv73B/3+c9/ntWrV3Po0CEA3n77bd55550BbcNAZcwki1WZYj7wHFAOrBaR14wx92P149gA/BD4kTHmTeB9rEAa4HPA3caYj4Fu4B9E5D033ojyhp1BBzpYzF12M2ZOLkM//fTTTJ8+nc7OTtasWcO6deuYOnVqbt6AcpUfs7zgzyoPXk0zrQpTKATr1weYNq2Fiy9uZfv2IOvXB6iuft32a3k109+NN97IV77yFR588EGeeuopbrnlFr761a/2DNx76qmnuOOOO/jzn/9MZ2cnd955J7W1taxZs4bZs2djjEk5cC9bK1asYNasWXznO9/pGbgH8IMf/IAZM2bQ1NTUZ+DexIkTef3113v+P8ceeyw//vGPOeWUUwa0HQNhrD7L/jF+/HjZulV7Y/iRNiSFz27/0dtuu41/+Zd/QUQwxvCVr3yFVatW5XQdKnt2/rebN28mHA737LtQKNSTyUn3+qXWz926lL6cNWsae973rFlLOOechXkvVan8ZdGi3vE1998Pr7/+Ouedd16+N0vZkGyfGWO2icj4ZMvrjHsqa34ZdKCcs5uRmzRpEmvWrKGzs5OKigomTZqUdnnt95w9twe8FUuW1211dbBsmXUpHaI9l9LXr8/3lik/STa+5n/9r3xvlXKbBskOFXr2xAm/DDpQ3pk6dSrr1q1j48aNTJo0KWNXC79NouIlvw1401q+2Ul1KV27hKm4VNWdfvObfG+ZcpsGyQ6UarYsPuggV4XQVWGYOnVq1v2QvZpExYvA2s2g18nJhNPMcCkcmwYqFII77giwZEmg5EpVqsxSja85coSe7kzK/5x0L9Yg2QEnDRwUR8bMq0EHqjDZzV46zTy73W/Wj7V8tbSZe7RUpUonWd/0UAh2767mwIEDDB8+XANlnxMRDhw4QHV1ta3naZDsgE78oFRqbk+i4qQ+9Nq1a3vWMXPmzLxXedCuEP7hdKIkpUaMGMG+fft49913870pKgvV1dWMGDHC1nM0SHZAJ35QKjecfJeGDBnSc9lMRBgyZEja5Xfs2EFXVxdgzSy1Y8eOnHdtKMUBb8VCS1UqpwYNGsQZZ5yR781QLtIg2SG7DZzdRtftAVDNzdao7sRGIBy2GgYte6TAu8+I3e/S4cOH097OhZqaGq666ip27tzJ6NGjNegtYqkupWuArAbKzrwCyp80SPaI2xM/2M0819X1vaSYeMlRKfDvZ2TUqFFUVFRk/d0YO3Ysr7zySs/yY8eOzbiO9vZ2/vVf/5Wuri727t3LKaecogGwUiprOq9AcdAg2UN2Mk1uD4CKX1K84452/v7v9/DLX45i3boazZ6oHvHPyPTpMG+eNaDJjX6advve2/1u1NTUMHPmTFe7Rzl5Hyr39AqZ8gudV6A4aJDsY24PgDr77HauvfZxRLq49tpyzj77ZkAbd9UrFLIC5PgsU9kEyHYuMTrte2+3a4Pb3aN0YK4/xK9+3HtvhIMHWxk6NMjixYG8X/1QpUfnFSgOGiQXCScDhzZt2sPevXvYu3c3I0eewaZNe5g5Uxt21ctuaSy7lxidllN0m9tXcnRgrjtCIStAvv32eoyJIlLJihUthEKawVPe0nkFioMGyUXETrYsHIb77z/A3r1r6erqpLy8ghdeuJqRI1MHQdqwlxYnpbHsXmJ0cgXEK25eySnlWutuO3iwNRYgd2FMlIMHWwENUJT3SnVegWI6TmmQXKLa2uCKK15nzZpuRASRbq644nXa2qamDID8mvVT7nBSGuu8886jrKyM7u5uysrKOO+889Kuo1gmyHB7YC7oSWq2hg4NIlLZk0keOjSY701SJabY+sa7Ofuo32mQXKIaGiASCfKTn/T2mZo9O0i6k16nWb9iOqssJU5KYw0fPpyZM2eye/duzjjjDIYPH55xPcVSOs3NgbmggwmzEQ7D4sUBVqxo6dMnubZWS7op7/i1MpATbs8+6ncaJJcwu32mnDTsxXZWqdIbNWpUz2fDb90n/MaPgwnjz3MrsHY7w9Z79SNAvItFba1ODKK85VVlICfsfr/dnn3U7zRILnF2+0zZbdiL7ayyFNk5qBZL9wk/cnswIbh/Ums3w2Y3qNaJQZRfOKkM5ITbXSG8mH3UzzRIVq5y2veyWL5ghc7JQbVYuk/4kdtlId3u0mE3w1ZMl61VabFbGQjst31edIVwEvQWUxuQVZBsjLkK+AFQDjwqIt/u93gV8DgwDjgA3CAie4wxVwLfBiqBKPANEdmUw+1XPmf3C6bdM/xFrwQULieNmxddOuxk2Px82VqpVOInc48+2s6JJ+5h3LhRTJ9ek/az297eztKlS3nrrbc466yzuOeee3J+tchpV4hiCnrtyhgkG2PKgYeAK4F9QJsxZoOI7ExY7FbgAxE52xhzI9AE3AC8B3xBRPYbY84HngNOz/WbUP5m5wumQZm77GYq3OxfVmwjwP3IyaQrbnfpCIdh/fp2Fi2yZvoMhdLP9OnVZWulMsn2+NnWZgXIr77aewL56KM309aW+rP+7LPPsnr1arq6unjhhRe48MILmTdvXtrtKfWuEF7IJpP8aeBNEdkFYIx5ApgCJAbJU4D7Yn8/Baw0xhgR+feEZV4DBhtjqkTkyIC3XBWlYuv07ydOu064dVDVS+n+5GaXjnAY7rijnenTrZk+p08v5447bubBB1MHD04uWyuViZtdGxoaYPPmvieQJ564h4aG5Ms3N8Mf/tB3+S1b9nDwYPqEQal3hfBCNkHy6UB7wu19wGdSLSMincaYPwPDsTLJcdcC25MFyMaYucBcgJEjR2a98ar4OA3KtB9zZk6z9G4dVOOX0qdNi3Dxxa1s3x5k/fqABkAFxO73ta0N7rprD3/8o/U5hC7uumtPygxb/MTprrueZteujdx11ySmT5+qXS7UgDhJGLjZtaGuDpYtm8qgQQ/S2RmloqKSX/1qKrNnZ34vGvS6y5OBe8aYWqwuGBOTPS4ijwCPAIwfP1682KZiFIlEimIKTLtfeu3HnB0/ZumrqyN89FE94bDVMFRXt6CzoxUWO9/XhgZobx/F44/3fg6vuGIUqZ5uBdVP09g4nc7OTioq1rBkybq0kx6BnjSXIjv73EnCwM2uDaEQrF8fYNq0TZow8JlsguS3gcS9OyJ2X7Jl9hljKoBhWAP4MMaMANYDN4vIWwPeYnWU5mYYMiRCQ0N9z8Qgzc0tHD4cKIm+naXcj9lP5dmc9DFevbqVzs4o0EVnZ5TVq9NPY60Kn53PYUMD3HbbRjo7OxEROjs72bVrI6tWTU35HD1pLnxuV3lwkjBwu2tDKAR33BFgyZKA9r33kbIslmkDzjHGnGGMqQRuBDb0W2YDMDP293XAJhERY8zxwP8H3C0iv8vRNhe95mYruEgUDlv3J1NXBwsWtHLkSJSuri6OHImyYEErdXXub6sfxA94xhhbZeY2b95Me3t7xmX9Kt4whMNhHn/88azeS01NDRMmTHAlaIj3MY5/duOXylN9DsNhePLJIFVVlZSXl1NVVcmTTwaP+uyr4mPnczhp0iQqKiowxlBRUcGkSZPSLp/spDkbxXBM8Cs7/1snxzW7+zwe8IZCIVsnUW4eP/v3vdfjoD9kzCTH+hjPx6pMUQ6sFpHXjDH3A1tFZAPwQ+BHxpg3gfexAmmA+cDZwCJjzKLYfRNF5J1cv5FiYndAUygETU1Bbr+9EmOidHdX0tQULJkz0VItM+e3DLrdcl1tbdYlxurq3lkfOzoCOjua6mPq1KmsW7eOjRs3MmnSJKZOnZp2eae12YvhmOAFP9bydZoZ9ss+Tmzj4wNTE2+r/MmqT7KI/Br4db/7FiX83QFcn+R5S4GlA9zGkuOkNuj8+QG2b29hzZpWZs0KMn9+aV2yLoYyc34qz+aUnXJdvV0w+s76qI2C6m/q1KkZg+M4J5fF3Z5EpVj4bcBbXKGXNuudTt26HY8BNGGQfzrjnk/ZrQ0aDsOvfhWgsTHQc6lGv1zJ+TG49Ft5Nqe0XJfyA7tZQi8mUSkGfhvw1v95hboPdDp1/9Ig2afsBBt6qcYeJwdht7NGfivP5oR+DlWh8mISFfBn9tnONvlxwJtSbtIg2YfsBht6qcY+OwdhL7JGfsxu26WfQ1XI3JxEBZwdR9wOqu1uUylmeVVp0yDZh+wGG3qpxl1eZI382HXCLv0cqlLhRb9nL07OnRzbNOD1l2KZH8GvNEj2IQ02/MVp1mjt2rU9z5k5c6Y2PkoVEbf7PTs9Od+2bRs7d+5k9OjRjBs3LqfbpPwlEolQX987P0JLS4sGyjmmQbJSGTjJGu3YsYOuri4Aurq62LFjhwbAOaTZk8JWivvP7nHESQC7bds2nn32WQB27doFkDZQLoYrWKWstbWVaNSaHyEajdLaqpMx5ZoGyUqpgqLZk8JWyvvPTvbZSQC7c+fOo25nyiarwhUMBqmsrOz5LgWDwXxvUtHRINnHSjHb4kdO+gaOHTuWV155pec5Y8eO9Whri59mTwqb7r/s2e3SMXr06J4Mcvx2On4cTKiyFwgEaGlp0TjBRRok+1QpZ1v8xunglpkzZ/qqzFyx0OxJYdP9lz27x4R41jjbPsl+HEyo7AkEAhobuEiDZJ/SbIt/OB3c4rcyc8VCsyeFTfdfdpweE8aNG5d1FwuvBhNqAkAVKg2SfUqzLe7yW3k2v06V7TfNzVBXB6FQb/YkHLbKIyarCtO7fO996ZZX3tDsV2ZeHBO8GEyoCQBVyDRI9inNtrjH6RTQbh7YnWarSy1DU1fXd2KdxIl3crG8Um5ye3Y7J9weTOgk2C+145ryLyMi+d6GPsaPHy9bt27N92aoIrZ582bC4TAigjGGUCjEhAkT8r1ZthuGUs3QxAPdefOsKdszTXttd3ml3FCqg+Tsvu9SPa6p/DHGbBOR8cke00yyKgp+zNDYZTdbXapdNEIhK+BdsgQaGzMHvHaXV8oNpTq7nd3ss/Z7zo52JfOGBsmq4NnNPBRLAX2/BvtuC4etjHBjo/U702yUdpdXyg2l+n0Fe8G+9nvOjtOuZFpa1h4NklXB0wxN6ZSZS2wI4sFu4u2BLq9Utux+l4rl5NxtXvV7LnShkHUcs9OVzG5pWc1Wa5CsioBmaEqnzFxbW9+GIN5QtLUlbxzsLq9UNpx+l4rh5NwLdv9PpdoG2O1KZre0rA58zjJINsZcBfwAKAceFZFv93u8CngcGAccAG4QkT3GmOHAU0Ad8JiIzM/lxqvipBkadxRDtiVZ9iJd94n48omXGEOhgAbIOVSKl2+L4btUTErxqhrY70pmt7Ssk2x1sckYJBtjyoGHgCuBfUCbMWaDiCROEn8r8IGInG2MuRFoAm4AOoBG4PzYT0krxcbELs3QuKdUy8zp7JXuKab/bTEM/i1lpXZVzUlXMielZUt94HM2meRPA2+KyC4AY8wTwBQgMUieAtwX+/spYKUxxojIR8BLxpizc7fJhamYGhM3aYbGPU6zLYXemOjsle4plv9tqQ7+LVXF0M447UpmdyKfUh/4nE2QfDrQnnB7H/CZVMuISKcx5s/AcOC9bDbCGDMXmAswcuTIbJ5ScIqlMXFCMzT+UYpl5nT2SvcUy/+2VAf/lqpiuKpmt+uZEzrw2ScD90TkEeARsCYTyfPmuKJYGhO7NENT2IrhpEVnr3RPsfxvi+FzrrJXqlfV7NKBz9kFyW8DiZ+EEbH7ki2zzxhTAQzDGsCnYoqlMbFLMzSFrVgGxNi9xKgy6y0P1fu/9Ut5KB38qzIpxatqdnmRrfa7bILkNuAcY8wZWMHwjcCMfstsAGYCEeA6YJP4bb5rHyjFhlozNIWv1AbEqOz4tTyUDv5Vbij0tkxrHjtTlmkBEekE5gPPAa8D60TkNWPM/caYa2KL/RAYbox5E/g6cHf8+caYPcD3gVuMMfuMMaNz/B58r7nZ+jAmCoet+wtNe3s7mzdvpr29PfPC9GZoQqGQBkwlIFm2RRUGu8epxPJQixb5p6+ifgaVG5y0ZXbbSzfFT2rj3/H4SW1dXX63y++y6pMsIr8Gft3vvkUJf3cA16d47qgBbF9R8GvGxS7N0KhMnE4pq5e588/Jccqr8lA6+Ff5QSFfVRtIzeNSLl/ri4F7xa5YCnKXYp8sZY/dvp1+a0hKmZPjlBfloXTwrypEfmwvnZzUlnr52ozdLVRfTrtOJH44583zT4Bs53JQPENjjNEMjUqppqaGCRMmZNUg6KVxf7FznErMNN9/f2+A3f/4OFBOPiN2PoNKucFJe+l294z+J7XZfFeTla8tJZpJtslp1wk/FuTWDI3KNz9fGi/FS4zhMDz4YIRQqJUHH0w/hbdX5aH8/BlRKhW/XVVzWvO4VMvXxhm/FaEYP368bN26Nd+bkVY4DNOmRbj44la2bw+yfn3qhiS+fKrAOp+B8ubNmwmHw4gIxhhCoRATJkzI3wapkuSkT7Lb/ZhL8RJj/LjW0VFPZ2eUiopKqqtbMh7fnLC7/7Tfuip2brfHA6luUewJA2PMNhEZn+wxzSQ7UF0d4aOP6gmHexsSSP3B8SrjYrch0QyN8gO7Azu96MdcijNktrXB9de3smaN9b4hyvXXt9LWltsg2cn+08G/qti5PQvgQGoel2L52jgNkh1YvbqVzs4o0EVnZ5TVq9M3oF4U5Hba8Gj3CVVovBgQU4qXGBsaIBIJ8pOf9L7v2bOD5Lpt9OOAJqXyTWcB9CcNkm0Kh+HJJ4NUVVX2XJJ88skgM2bkt+uE04ZHMzSq0HhRZq5UZ8h08r71CpZSuaGzAPqPBsk2tbXB+vUBqqt7G5KOjkDe5zLXhkeVCq8GxJTqJUY771uvYCmVP9ruu08H7vmYDm5RauB0gKp79H+rVH5pnDBwOnCvAOngFqVyw262xe4o8IGMGi90mslSKr/8NAtgMR4LdTIRn9JJFpTKjfjl/VAolFWjEK+F/swzVmH/Z55pZ/p06/50y8cL88dLPKZa3gtOJz0CexMa2P3fKqXyx+24wo/HwoHSTLKH7Fzm0AyNUrljJ9sSCsGjj7bT1vY4FRVddHaW8+ijNxMKJX++H6eddzrpkV7BUqp4uV1mLn4stDOPhN9pkOwRnd1OqcJx4ol7qKjowhihoqKLE0/cA6RvHOLTOTc25n82TaeBu46WV6p4eVFmzu48En6n3S084uQyR01NDRMmTNBGSimPvf/+KDo7yxExdHaW8/77o9IuHw7D+vXtLFq0mfXr24/q6pAPVuDeTji8mXnz2rMK3OOZJmOMXsFSqgjZjSvsxi7J5pEoZBoke0QbH6UKQzgMc+bUUFd3M/X1IerqbmbOnJqUgW84DHfc0c4llyzld79bzCWXLOWOO3IbKDvpY/zMM+3s3buUysrF7N27lGee0T7GSil77MQuifNIlJeXU1VlzSPhh6SBU9rdwiG7ZVS0+4RShaF3Gvka4l0sjjsu9TTybW0wdeqzfPvbq2OXJF/g7rsvpK1tXs66XdjtYxwOw/Llz7JtW+82/ed/Xshxx2XeJu1jrJSKsxO7xOeR+POff8rGjRuZNGkSw4blfx6JgdAg2QGnZVS08VHK/+xOI9/QAAsW9L0kGY3uyWnJI7t9jNva4Jxz9vDyy73bdM45ewq6sVJK5Ue2sUtDQzw+eo3TTz+d1157jZtvHpdy0HMhyKq7hTHmKmPMG8aYN40xdyd5vMoY8/PY4783xoxKeGxh7P43jDGfz+G259bzz0NFBWzalHFRLc+mlEo0depUqqqqKCsro6qqiqlTp2Z8TtpSa0mOR3b6GDc0wD/8Q99t+od/mFqwtUqVUoXBUXxkI/7yWsZMsjGmHHgIuBLYB7QZYzaIyM6ExW4FPhCRs40xNwJNwA3GmNHAjUAtcBrwW2PMJ0WkK9dvZMBuuAG6uuC66+D999MuquXZlFKJAoEAmzZt6pmqPtO0zhmvRiU5Hj3zTDvd3Y9TX2+VpXvmmZuZMiV1hsbuNiml1EA5io9sxF9eyzgttTEmANwnIp+P3V4IICLLE5Z5LrZMxBhTAfw/4GTg7sRlE5dLtT7Pp6U2JvVjaf43OrWjUsqplNM5pzgeCTDp8y8SCIRjtwyRSIgFCyZo9wmllK9kHR85jL9yLd201Nl0tzgdSLweuC92X9JlRKQT+DMwPMvnYoyZa4zZaozZ+u6772axSTn0m9/AkCF97zvmGGhpSfs0Lc+mlHIq5YjxFMejx76yhXnzRlFRYT2noqKcefNG0dbm+aYrpVRaWcdHDuMvL/li4J6IPAI8AlYm2dOVX3klVFXB4cO991VWwhVXeLoZSqnSkXLEeIrj0axHrK4SF1+sFXKUUkWiAOKvbDLJb9N3qqkRsfuSLhPrbjEMOJDlc/Pv8GE44QRoarJ+J+4wpZTKoXjN48RsS5+ax2mOR3oFSylVVHwef2UTJLcB5xhjzjDGVGINxNvQb5kNwMzY39cBm8Tq7LwBuDFW/eIM4Bzg5dxseg51dFidxRsarN8dHfneIqVUkYrXPI4X2I/XPK6riy2gxyOlVKnw+fEuY3cLEek0xswHngPKgdUi8pox5n5gq4hsAH4I/MgY8ybwPlYgTWy5dcBOoBP4mi8rWyillEfs1jxWSimVHxmrW3jN8+oWSimVB7NnR1izppVZs4KsXp1debZIJKIl3ZRSKofSVbfwxcA9pZQqJStXRlizph5joqxZU8nFF7cwf376oDcSiVBfX080GqWyspKWlhYNlJVSykVZzbinlFIqN8JhWLCglbKyKCJdlJVFWbCgtaePciqtra1Eo9HYtNdRWltbPdlepZQqVRokK6WUh9raoKkpSFVVJeXl5VRVVdLUFMxY8zgYDFJZaT2nsrKSYDDoyfYqpVSp0j7JSimVB076F2ufZKWUyq10fZI1SFZKKaWUUiVpoNNSK6WUUkopVVI0SFZKKaWUUqofDZKVUkoppZTqR4NkpZRSSiml+tEgWSmllFJKqX40SFZKKaWUUqof35WAM8a8C/wxT6s/CXgvT+tW3tP9XVp0f5cW3d+lR/d5acnV/v6EiJyc7AHfBcn5ZIzZmqpWnio+ur9Li+7v0qL7u/ToPi8tXuxv7W6hlFJKKaVUPxokK6WUUkop1Y8GyX09ku8NUJ7S/V1adH+XFt3fpUf3eWlxfX9rn2SllFJKKaX60UyyUkoppZRS/WiQrJRSSimlVD8aJAPGmKuMMW8YY940xtyd7+1RuWeMWW2MeccY82rCfScaY543xvx37PcJ+dxGlTvGmBpjTNgYs9MY85ox5h9j9+s+L0LGmGpjzMvGmB2x/b04dv8Zxpjfx47tPzfGVOZ7W1XuGGPKjTH/box5NnZb93eRMsbsMcb8wRjzijFma+w+14/nJR8kG2PKgYeAScBo4IvGmNH53SrlgseAq/rddzfQIiLnAC2x26o4dAL/JCKjgUuAr8W+17rPi9MR4AoRGQtcCFxljLkEaAL+j4icDXwA3Jq/TVQu+Efg9YTbur+LW0hELkyojez68bzkg2Tg08CbIrJLRKLAE8CUPG+TyjEReRF4v9/dU4C1sb/XAlO93CblHhH5HxHZHvv7IFZDejq6z4uSWA7Fbg6K/QhwBfBU7H7d30XEGDMCuBp4NHbboPu71Lh+PNcg2Wo42xNu74vdp4rf34jI/8T+/n/A3+RzY5Q7jDGjgIuA36P7vGjFLr2/ArwDPA+8BXwoIp2xRfTYXlweABqA7tjt4ej+LmYC/MYYs80YMzd2n+vH84pcv6BShUhExBij9RCLjDHmWOAXwJ0i8hcr2WTRfV5cRKQLuNAYczywHjg3v1uk3GKMmQy8IyLbjDHBPG+O8sbnRORtY8wpwPPGmP9MfNCt47lmkuFtoCbh9ojYfar4/ckYcypA7Pc7ed4elUPGmEFYAfJPROSXsbt1nxc5EfkQCAMB4HhjTDwZpMf24nEpcI0xZg9WF8krgB+g+7toicjbsd/vYJ0EfxoPjucaJEMbcE5sVGwlcCOwIc/bpLyxAZgZ+3sm8Ewet0XlUKx/4g+B10Xk+wkP6T4vQsaYk2MZZIwxg4Ersfqhh4HrYovp/i4SIrJQREaIyCisNnuTiNyE7u+iZIw5xhgzNP43MBF4FQ+O5zrjHmCM+Tus/k3lwGoR+VZ+t0jlmjHmZ0AQOAn4E3Av8DSwDhgJ/BGYLiL9B/epAmSM+RywGfgDvX0W/zdWv2Td50XGGHMB1sCdcqzkzzoRud8YcyZWpvFE4N+BL4nIkfxtqcq1WHeLu0Rksu7v4hTbr+tjNyuAn4rIt4wxw3H5eK5BslJKKaWUUv1odwullFJKKaX60SBZKaWUUkqpfjRIVkoppZRSqh8NkpVSSimllOpHg2SllFJKKaX60SBZKaWUUkqpfjRIVkoppZRSqp//HxdMWFa3ABWzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "times = mod.two_layers_forward_refraction(sources, receivers)\n",
    "plt.subplots(figsize = (12, 4))\n",
    "mod.plot_response(sources, receivers, obs_times, marker='x', color='blue', label='picked times')\n",
    "mod.plot_response(sources, receivers, starting_times, marker= '.', color= 'grey', label='starting model')\n",
    "mod.plot_response(sources, receivers, calc_times, marker= '.', color= 'black', label='fitted model')\n",
    "# plot legend\n",
    "handles, labels = plt.gca().get_legend_handles_labels()\n",
    "by_label = OrderedDict(zip(labels, handles))\n",
    "plt.legend(by_label.values(), by_label.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "nbsphinx-thumbnail": {
     "tooltip": "This tooltip message was defined in cell metadata"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-2-13a568e85554>:68: UserWarning: linestyle is redundantly defined by the 'linestyle' keyword argument and the fmt string \"-k\" (-> linestyle='-'). The keyword argument will take precedence.\n",
      "  plt.plot(x_span, z[i], '-k', *args, **kwargs)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7faa2fdff6d0>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsEAAAD4CAYAAAANWzs4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAABDm0lEQVR4nO3dd3zb1b3/8df5atrWiGecxM5OSMIuYUNZYbeMQrmlg9KW0t7bvaF0UCgt3XT30j24pdwO6A9KoXChLdACYRVoGElIJMd24ji25W1LOr8/NCzZkjPsxI71fj4eedj+fo+OjvX9Rvrk5JzPx1hrEREREREpJc5UD0BEREREZF9TECwiIiIiJUdBsIiIiIiUHAXBIiIiIlJyFASLiIiISMlxT8WT1tTU2IULF07FU4uIiIhICXniiSe2W2trRx+fkiB44cKFrF27diqeWkRERERKiDFmc6HjWg4hIiIiIiVHQbCIiIiIlBwFwSIiIiJSchQEi4iIiEjJURAsIiIiIiVnUoJgY8xZxpgXjTHrjTFXTUafU66lBU46CVpbx212+OGHY4wZ8+fwww9Xe7Uft/10GovaT0H7Ue8xUz6efdh+Oo1F7dV+Iu2n01j2RfusXYyRprsJB8HGGBfwXeBsYBVwqTFm1UT7nXLXXw8PPQTXXTdus2OPPRav15t3zOv1ctxxx6m92o/bfjqNRe2noP2o95gpH88+bD+dxqL2aj+R9tNpLPuifdYuxkjTnbHWTqwDY44FrrXWnpn++WoAa+0Xiz1m9erVdl/nCf7gBz/Iz372szHHa2trmTdvHolEgmeffZZDurry/mVwefrPFmM4MBQq+PiNGzeSTCbzjodCIYwxNDQ0UFNTQ19fHy+99BLJZJLu7u68to7jcNBBB7F589g0dg0NDaxbt65o/0uXLiUQCNDR0cHmzZuL9r9s2TJaC/yLbfHixTzzzDNF+z/wwAPxeDy0trbS2tpatP8FCxawY8eOMf2vXLmSxx57rGD/LpeLQw45BIDNmzfT0dFRtP/6+np6e3vzjvt8PhYvXly0/4qKCpYvXw7ASy+9RF9fX9H+q6qqGB4eHtPHvHnzivZfVVXFggULAPjXv/5FMpks2n8gEMAYk3e8srKSSCRS9LWvr6+nvr6e4eFhnn/++aJ9H3bYYWzYsIHR6uvrefnll4v2v2DBAiorK+np6WH9+vVF+1+1ahXRaHRM//Pnz8+Oq1D/y5cvp7y8nO3bt9PU1FS0/yVLlrBt27Yx/S9dupSnnnqqaP8HH3wwLpeLLVu20NbWVrT/xsZGOjs7xxw/4IADil5br9fLgQceCMDGjRuJxWJF+6+rq6O/vz/veHl5OQsWLCjY/9HAMcBN6Z/fDDQBA8Cjo14Dx3GYNWsWiUQi73hlZSX19fVFx19TU0NjYyMATz/9NEDB8RtjCAaDY+7N+vp6KisrefTRRxn9GZH5e1FXV8fAwAAvvPDCbvc/d+5cXnzxxaLXdtGiRYTDYbq6unjllVeKvvYHHHAAzc3NjLZw4UKeffbZov2vWLECv9/Ptm3baG5uLtr/okWL2L59+5j+ly9fzhNPPFG0/8MOOwyAaDRKe3t70f7nzZtHLBbLO+7xeFi6dGnRa+v3+1m5ciUA69evp6enp2j/NTU1DA4O5h0PBAI0NjYW7T8cDrN48WIAnn/+eYaHh4v2HwqFxtwf1dXV1NXVFe2/rq4u7zMXdv/eCYfDRe/NxsbGvM/c3e1/0aJF+P3+ov0vW7aspD9zI9u340sm+SDwdM65JLB+zpwxn7m5Lr/8cm666aYxY94XjDFPWGtXjz4+GcUy5gG5n5BNpN7nRw/gSuBKSH14ToWheJJ4Mv+mbu0aoMvEsMkkfUMJnvCXs3BoiOpkHAcYBP43UMmnq2bT19Y0ps+tPXG8gUoGYu3ZY8blpn84dQNGd/SzbShGfGiAvqFE9rxNxLPtvcFKop2D2fO5tnQNjdv/xrYe3LEkQ3294/bfEhumv0D/mzvG7//F1m4cl5uBWD8D4/Tf1jPMYIH+N7QPFu8/bvl3c+oDoDc2yHDB/g3eYCUd/Qnio/ofSA6P2/9Q71C2/+7eIRLDxfvvHkyQjOf3P9QzxMA4/W+LDdKb7r9nMA7pN8xC/fcP5X8QASR6x793mjv72ZGMkUzEx722m3YUvndaYsPj9r+5vZeWfhfxwb5x+28qcm827eTeXL+tB7c3zmBPf/beK9R/a/dw9t7KtWkn9+YLLd0Yx6G/ayB77xV67bf3DjM0un+THPfeGUgMZ++dnp6h7L1XqP+ugbH35qAdIl6k/2c9XqoTCfriQ5RbSxxoc9xs8noxw0Nj+u8ZjGOT+f0Pxwbp9RQf/9bYAN3p8edeu4L35vDQmNc+de+V462YxWBPx5j+mzr62R6PkRgeGrd/T0WY/uE4ozXHxr+2m7b34uk1DPeP/77WHBsqeG9GOsfv/+WtPbg8Qwx0j/++tq2n8L35yo7x+8/cO32xgey9V+i1b++NZ9/3sv3Ex3/fHLQ592bv+PdmbCCRfd/LGOodYmic/od7hhhI9x/rH87ee4X67x0ayr7vZcRjA3S7ivc/+jM393xe/4FZ9I+amADS954fb0WYwZ7OMf0X+swt1L+nIkT/cIH3ne29eMoMnrIgQ32xvMf3DydL/jP3iAWruLZtC8Pp94Uk0J5+/0oW+czNPr7A6z3V9lnFOGvtzcDNkJoJ3lfPm3HTTTfhPeFt3LZ2bCCbUZH+evU93+WNT/+ZIZcbbyLOLUuPYeDM9zC3yOPiPTto/u8rsPEhjNvLvHf9GFegcky7UJH2NW++CVegMnt+zLh2of9yYNZO+h87opSyXey/2O+b6b8Y74T698zo/nf13gkUaZ/pO1hk7JN1b4aL9F++i69N5sju3jv+/fjajtf/pfd8F//Tf2bA5eF/EnFuOeR0Pn3me/b5+Gvf8s0J9597743uv+6t357wvR8u0j7z2swqMvZdvTeLjWdnr71vht6b+7L/ipzvx9ybl31rwv3nfqaOvTe/M6H+S/UztwcYuOe7fDs3Rkq/fxXrP+MzV51adMxTZTI2xm0BGnN+bkgf22/V9HXyq8PP5sLLvsavDj+b2r7Ocdu7A1VUHHQaYKg4eM24N6faq/10HYva77v2xd5j9pfxT0b76TQWtVf7ibSfTmPZF+13N0aaziZjJvhxYJkxZhGp4PcNwBsnod8p8+4Lr8l+/5kz/muXHhM+/lKG26PMOu5StVf73Wo/ncai9vum/XjvMfvD+Cer/XQai9qr/UTaT6ex7O32exIjTVcT3hgHYIw5h9Q+DxfwE2vtDeO1n4qNcQAf/+0z4y6HEBEREZHJ98hVpzJ3VtmUPPfe3BiHtfZPwJ8moy8RERERkb1NFeNEREREpOQoCBYRERGRkqMgWERERERKjoJgERERESk5CoJFREREpOQoCBYRERGRkqMgWERERERKjoJgERERESk5CoJFREREpOQoCBYRERGRkqMgWERERERKjoJgERERESk5CoJFREREpOQoCBYRERGRkqMgWERERERKjoJgERERESk5CoJFREREpOQoCBYRERGRkjOhINgY83pjzPPGmKQxZvVkDUpEREREZG+a6Ezwc8DrgL9NwlhERERERPYJ90QebK1dB2CMmZzRiIiIiIjsA/tsTbAx5kpjzFpjzNq2trZ99bQiIiIiImPsdCbYGHMfUF/g1DXW2jt29YmstTcDNwOsXr3a7vIIRUREREQm2U6DYGvtmn0xEBERERGRfUUp0kRERESk5Ew0RdqFxpgm4FjgLmPMPZMzLBERERGRvWei2SH+APxhksYiIiIiIrJPaDmEiIiIiJQcBcEiIiIiUnIUBIuIiIhIyVEQLCIiIiIlR0GwiIiIiJQcBcEiIiIiUnIUBIuIiIhIyVEQLCIiIiIlR0GwiIiIiJQcBcEiIiIiUnIUBIuIiIhIyVEQLCIiIiIlR0GwiIiIiJQcBcEiIiIiUnIUBIuIiIhIyVEQLCIiIiIlR0GwiIiIiJQcBcEiIiIiUnImFAQbY75ijHnBGPMvY8wfjDGzJmlcIiIiIiJ7zURngv8CHGStPQR4Cbh64kMSEREREdm7JhQEW2vvtdbG0z/+E2iY+JBERERERPauyVwT/Hbg7mInjTFXGmPWGmPWtrW1TeLTioiIiIjsHvfOGhhj7gPqC5y6xlp7R7rNNUAcuKVYP9bam4GbAVavXm33aLQiIiIiIpNgp0GwtXbNeOeNMZcDrwFOs9YquBURERGRaW+nQfB4jDFnAR8HTrLW9k3OkERERERE9q6Jrgn+DhAE/mKMedoY84NJGJOIiIiIyF41oZlga+3SyRqIiIiIiMi+oopxIiIiIlJyFASLiIiISMlRECwiIiIiJUdBsIiIiIiUHAXBIiIiIlJyFASLiIiISMlRECwiIiIiJUdBsIiIiIiUHAXBIiIiIlJyFASLiIiISMlRECwiIiIiJUdBsIiIiIiUHAXBIiIiIlJyFASLiIiISMlRECwiIiIiJUdBsIiIiIiUHAXBIiIiIlJyFASLiIiISMmZUBBsjLneGPMvY8zTxph7jTFzJ2tgIiIiIiJ7y0Rngr9irT3EWnsYcCfwmYkPSURERERk75pQEGytjeX8WAHYiQ1HRERERGTvc0+0A2PMDcBlQBdwyoRHJCIiIiKyl+10JtgYc58x5rkCf84HsNZeY61tBG4B3jtOP1caY9YaY9a2tbVN3m8gIiIiIrKbdjoTbK1ds4t93QL8CfhskX5uBm4GWL16tZZNiIiIiMiUmWh2iGU5P54PvDCx4YiIiIiI7H0TXRN8ozHmACAJbAbePfEhiYiIiIjsXRMKgq21F03WQERERERE9hVVjBMRERGRkqMgWERERERKjoJgERERESk5CoJFREREpOQoCBYRERGRkqMgWERERERKzkTzBO83HnjgAX57w6dptwFcoVrc6T+uUB3uYDXG5ZnqIYqIiIjIPlIyQXBfXx+VPevpaW2nrSeed84Ygy8QwgTrMMG6/CA5mPrqlIcxxkzR6EVERERkMpVMEHzuuedybvwieOpX9A9bmmJJIl2WSFfmax/R2EYiOzYR2ZigeziZ93i3x4M3VEMykA6Sg7WjZpRrcDz+KfrtRERERGR3lEwQnKvMY1hW7WJZdeHz1lra+21OkJwk2mWJxNqIdLUR2Qwt3XGszX+cryKIK1iDDc7Om0XOBMmuikqM49r7v6CIiIiIjKskg+CdMcZQU26oKYdXzSkctA4lLFtilmgsmTObPECkK0ok1kQkkqBzMH822XEc/OFqbKAOJ5SZUa5JB8m1uEN1OL7yffErioiIiJS0kgqCb3toPc89NkBjyGF+2KExbGgMOQR9u7/W1+syLKo0LKosnmCja2BkJjkzqxyNdRLp2kGk+SWa1sVJJPOnk73+MtyhWmywLjuTPLLsog5XoArjKqnLJiIiIjLpSiqauv9fTfzwb0Pkhp2VftjxiRAAX3pokI0dSRrDDo0hQ2PYYeEsh8XjBLrjCfsNB/tdHDy78GxyImlp6bFEu3LXJ8eJxJqJxJqJtkJ779hNfP7gLEhv4ssPklPfO/6gNvGJiIiIjKOkguD//q+T+c5RUZq7U8sYol2W/vhISPxCe5I7X4qzvW/k2CGzHZ55dwCAd9zRT3u/zQbIjSGHFTUOhxdZMrEzLsfQEDI0hByObSzcpncos+Qid1a5h2gsRqRtI5H1cYbi+bPJbq8Xb6iWZCA/SM4Gy8EajNu7R2MWERERmQlKKggG8LgMC2YZFswaO7v70/PLALLZI6Ixm7f5LQls6Ejy4KYkXYOpY+cuc3PnG1PreI/+UQ8Gk11mMT9sOHKui+Pnp17mpLU4uzlDW+E1rKhxsaKm8PmktbT12pzlFpmAeSuR2FYir8DW7viYx/kDYUywNi8l3EiQXItTEcYY1VIRERGRmankguBdUSx7RCZIBugeTM3Q5jp0touNHUn+tTXJXS/F6Y/Du47wcPx8N4mkJfjFbmYHUgFyJlA+a6mbkxe6SVpLR7+lqszs1lIGxxhmBwyzA3DkvMIz0oNxS1Ns9PrkPiKxV4h2bmbzpgTdQ/m/i8vlxheuTqeES88opzfxuUN1qZRw3rKCzyciIiIy3SkI3kNBn2FVbX7QefNrR4JCay07+i2J9EzyYAI+cLSXaDoYfSSaYEssTtBrOHmhm5ZuS8M3eihzk7cm+fJDPZy00E3/sM2uVw7t5kY+n9uwpMqwpKrwzK61lo4B8oLkaFeSSKydSFc7keg6mmNxkqNTwpVX4ArVYQM5M8nB2pGgOaCUcCIiIjI9KQjeS4wxVJePBKvlHsMX1+QX00hay3Ai9X2Zx3DTmb5UAJpeivGXDXHWLEpdoqdbExz3kz4AQj6ys8mfPcnHMQ1utvUmeXZrMjvDXObZ9UDZGENVGVSVuTisvnDQOpyw2bXUI+uTh1Ip4bqbiG5J0tmfyO/XcfCHqiBQh8nNcJGdUa7F+Cq0iU9ERET2uUkJgo0xHwG+CtRaa7dPRp+lwDEGX/oKVJUZPnCMr2jbpVUOv76ojGhXKkBObexLZmdn/7opwSW/7c+2ry5LrU3+xQVlHDzbxQvbEzzZkszOMM8LGjyuXQ8+x1tLnREbzGS6yE0J10Wkq4NI68tEX4gTHzWd7PH58YRqSQbrcI9KCZeaWa7GuDy7PE4RERGRXTHhINgY0wicAUQmPhwpprbC4Q0HFQ9AT1vs5v8uK89mvcjMJmeWTvzp5TgfuXcw294A9QHDo1dU0Bh2+OumOE+0JFL5k9OBcn3A7NZGvpDPcGCdiwPriqeE29prRyrwZQLmWAuRWAvRbdDWUyAlXCBcMCVc5qtTFtJssoiIiOyWyZgJ/gbwceCOSehL9lBVmeGURcUv57tXezl7qTu7Jjkzo1yTXrJx50txvvqPobzHeBzovjqIz2341b+G0sstnJwUcYbail3PIOFyDHODhrlBh2MaCrfpS2fmyE8J10s0toFI+ytENiYYGB61ic/twReuIRmswxWswx3KVOEb2dDneIrPsouIiEjpmVAQbIw5H9hirX1GM3HTW7nHsLLWxcrawue/fLqPq0/0pYPj1ExtW5/F505d10eiCX705DC58WdNuaHtY0EArvtrutBIToC8qNJhRc3ubYwr9xiWV7tYXl34vLWW7X2FUsK1EYm1Edn8PK3d8bzUdgD+iiBOzia+0UGyK1CplHAiIiIlZKdBsDHmPqC+wKlrgE+SWgqxU8aYK4ErAebPn78bQ5R9IXdz3KEFNsd979wyvnOOn229Nrvcon945Hy0K8n9r8Rp7rbZdcqH1zs8+a5UoZE3/76f9v5kakNfelPfqloXRxVJ6zbeOGsrDLUVcMTcwo8dSli2jEkJN0AktoloLMLmSIKOwfzZZMflwh/KpITLBMg16WwX6WUXvvLdGquIiIhMXzsNgq21awodN8YcDCwCMrPADcCTxpijrLWtBfq5GbgZYPXq1Xb0eZn+HGOoDxjqC+Qk/uF5qfRw8WQ6i0TOpj2Acg+8sN3yZEucbb2pE+cd4OaON6QCy0O+34PLIVtkpDHkcGyji1cvcGf7dTu79r8NXpdhUWVqJroQay1dg4VSwu0g0rWDyJYX2LIuQWLUJj5vWTnuUC02MDKLnLdGOVCtlHAiIiL7iT1eDmGtfRaoy/xsjNkErFZ2iNLmdgzzw4b54fwANDeH8kA8NVObTK9ZsNZy4nwXm7osmzqT/D2SpHMA/mu1h1cvcDOcsJR/oZva8pFqfI0hh9ce4ObURalCI609ltkVBtcuBMrGGGb5YZbfxSGzCwet8WSqz9xMF5GuYSJdW4jEthBtsezoS4zp1x+qzNvEl58/uRbHH9AmPhERkWlAeYJln/Oni3dkGGP47rn51ed6hixD6RhzOAnXnOjNLsN4vi3Jn9fHmR0wnLrITVPMsuCmHtwOzAuOrEm+8ggvJy900ztkebE9tV65pnzXKvK5HUNDyNAQcjiusXCbnqFMSriRpRfRWDeRri4i2zYQfTnOUHxUSjivD0+4hmRgdk4VvpysF8EajFsp4URERPa2SQuCrbULJ6svkYA3v9DItSfnFxqxdqQaX9Br+N45/mxauGhXkke3JLhgRarB060JTvhpqtCI3w0N6SUXnz/Fx7GNblp7kjzVkkgHzw5h/67N1Aa84282TFrLtl47stwiEzDHtqb+bIBto1LCAfiDYUygDpNdj5ybGq4Gp3yWZpNFREQmSDPBsl8yxpBOXEFlmeE/j/QWbbuixuF3l4wuNGKzSyce3JTg0t+NFBoJelOlq2+9KFVo5PltCR7bkshLD1e+CxX5ctdQF9sAOBC3NI3ZxNdHNLaRSMcmIpsSdA+NTgnnxhfKTQlXiys0UoXPFazF8foLPp+IiIikKAiWGa+63OF1K4unPzt7qZtH3l4+UrI6/bWybKTQyMfvG8x7THWZ4el3V9AQcrhvY5y1zYls1ovGkMO8kMG7CxX5/G7D0irD0qrim/h29OemhMsEzNuJxLYTif6b5tjYlHC+8gCuUC02kF9gJPO9q6JSm/hERKSkKQiWkhf2G45tdHNskbW/7z/ay0WrPHlFRqJdSWrThUbu3RDnK4/kFxpxDPRfE8TrMvz0qSGe2ZqfQ3l+2GFeaOd5iY0xVJcbqsvh8DmFg9bhRCojR/4mvkEisSjR7iYiW5J09udv4nMcB1+oCpueTR5dhS+VEq5iF149ERGR/ZOCYJGd8LkNiysNi4ukXPvy6X4+c5KPppxZ5O19NjsT/GRLgp8+PUxvTl7l2RWG1o+mCo186v8G2NCRn0N5SaXDwUUyV4zmcRkWzDIsmFU8qO4asOnCIrkp4bqIdHUQaXmZphfixEenhPOXpVPC1aWXW9TlzygHqjEuvYWIiMj+SZ9gIpMg4DWsqHGxombsuW+fU8a3zvbTOUB6uUWS/pz9cO19lse3JPj9ung2I8bquQ6PvzNVaOSi2/rY0W+Zn1mTHHI4eLbDcY27/tc37DeE/S4OqiscWCfSKeGieSWrE0S6monEmoluhe29+Zv4jDH4g7PyUsJlZ5SDqTXKTllIm/hERGRaUhAssg8YY6gsg8qysbmJv/+aVHq4pLW09VqiMZs3Kzu7wtDaY3kgXZEvYeHCFW5+/x+pv77Lv51KD5ebQ/mE+S5OW5w63z9sKdvJRj6XY5gXMswLORzTULhN33ChlHA9RLpiRLZvJLI+zuColHBujxdvuIZkoC41k5wOjnNnlI27+KZGERGRvUVBsMg04RjD7IBhdiD/+PdycihninjE0wkjrLWcu8zN5q4k0ViSZ1rjbO21vO8oL6ctdjMYt1R8oZvKMpO3JvnCFR5OX+ImkbRs7rI07MJGvnKP4YAaFwcUmO3OjKWtr1BKuG1EurYReQVauwukhAuEMMFaCNbhDo5em1yHUxHGmJ2vnxYREdkdCoJF9iOZIh4Zxhi+cVZ+OrTBuGUwvawinoQbTvVlcyhHupI8HEmyaJbD6UvcRLosS77VA0B9YCRQfs+RXk5d5KZ70PLctlR6uDmB8SvyGWOoqzDUVcDquYWXXQzGLVu6R6eE6yca20SkM8LmTQk6RqWEc1wu/KFqksHZedX3coNlx1tW8PlERESKURAsMsP43AZf+m92hddw9Ym+MW0yJatn+Q0/Ps+f3dAXjSVZ15akcyB1/qnWBCf9LFVoxGVgbroi31dO93Fco5stsSSP56SHq91JRb6dbTK01tI5QM5yi0zAvINIbAeRpnVs6YqTHJ0SrqwilRIuWIcrmJ8OLrWJr0op4UREJI+CYJES5KQD1coyw9sPL74m96A6F3deWpZNC5cpNuJzjRQaefMfRgqN+FzQEDL84T/KOXi2i2daE/yjKZGdYZ4fdgj7KBoo566dPrS+cNAaT1paxqSEGyISayLavYVIc5KOvvyUcMZx8Acr8zbx5Vbhc4XqcHwV2sQnIlJCFASLSFFVZYZzl3uKnj/vADdr31mRV2Qk0pWkJp1D+e71ca6+P7/QSMALL743wNygw90vD/PYluTIpr701wpv8WDU7ZhUu7DD8UXadA/mZ7pIpYSLEenqJLJ1PdGX4gwn8qeTPT4/nlDOJr7cKnyhOtzBaoyr+GshIiL7FwXBIrLHgj7DEXNdHFFkDfDHjvPy5kM8eUVGojGbLTTy4KYEXy5QaGTgmiAel+G/1w7xzNb8anzzww6LiiynyB3XqloXq2oLn09ay9Yem7PcIvOnlUislegG2NZTICVcIAzBWkxwVM7kUB3uYA1OeVizySIi+wkFwSKy17jSG/kaQg7HFjj/pdP9XHeKjy3dmQA5SXufxZNebrFue5LfPB9nR//IrO3coGHLh1OFRj5270ihkfnp2eGlVQ6vKlJdL8MxhjlBw5wgHDWvcNv+YUtTXt5kSzTWS6Srm8iOTUQ2Jugezt/E53Z7sinhRnIm5wbLNTgef8HnExGRfUtBsIhMqfE2y910lp+bzvLTO5QKSKMxS//wSEA8ELe82J7kLxvj9KQnlI+e5+KfV6RKPr/mf/roGLDZIiONYcOhs12ctDD11metLTpzW+YxLKt2say68LittbT356aESwfLsTYiXW1ENkNLdxw7ehNfRRBXMLWJb2yQXIsrUKmUcCIi+4CCYBGZ9iq8hXMUf/ucVGo0ay2xwVRFvuGcPXFLKh2e3ZbgiZYkt78QZzABF610Z4Pg+Tf14HWRDpBTOZRPXujmjCWp850DtuhGPmMMNeWGmnKKzjwPJSzNY1LCDRCNRYjEomyOJOkczN/E5zgO/nA1NlCHk1uqOliT/r4Ox1e+py+liIikKQgWkf2eMYawH8L+/GD0m2ePLD2w1rK9z2ZLU1tredPBnmxw+rfNcbbELP3DcMYSN/3DlsovdVPhIRsgN4YcLjnQw5lL3cSTlpfbkzSGHQJFNvJ5XYaFswwLZxWf2e0asAVSwnUS6eog0vwSTeviJEblhPP6y3CPSgmXW2DEFajCuPT2LiIyHr1LikhJMMZQW5FfaOTGNfnrcxPJkUIjSQtfO8OXl0P5ufVxDp7tcOZSN5s6Lau+1wtApT8TKDt86JhUtb6uActTralNfQ0hg89dOFAO+w0H+10cPLvwbHIiXSUwPyVcnEismWismUgrtPcW2MQXnJWXEi4zi5xdm+wPahOfiJQ0BcEiImkux1CenrSt8Bo+fOzYQiMZNeWGX11Ylk4PN5JDuS+9ZvnJlgSn/qIv276uIlWR71tn+zmu0U2kK8k/ools1os5QYO7QEU+l2OYFzLMCzkc21h4LL1DhVLC9RDpihFp20h0fZzBeP5sstvrxRuqTaeEG1uFzx2swbiL55AWEdnfKQgWEdkDs/yGNx1SPG/w4XNc3PeW8mwO5Ug6+0WFJ5MeLs5bbx/ItncZmBM0/OmNqUIjTzQneCgSzxYZaQylZrKdArO3FV7DihoXK2rGnAJSKeHaekenhLNEurYSiW0l+gq0dsfHPM4fCGNCtZjcbBc5qeFSKeG0iU9E9k8TCoKNMdcC7wTa0oc+aa3900QHJSKyv5vlN5y2uPhb7MWrPLxqjis7i5xZE1yXXrJx74Y4n/y//EIjXhds+kCAOUGHP744zD+b8nMoN4YdKv1jN/I5xjA7YJgdgNVFcjoPxi1No4LkaKyPSNcrRDo3E9mUoGMoPyWcy+XGF65Op4TLmVEO1mQLjjjesj15+URE9rrJmAn+hrX2q5PQj4hIySj3GA6qc3FQXeGg9KoTvFzxKk9+yequZDZI/kc0wVceGSK38J3LwOCngrgMfPvRIZ5uTeQFyAvCqSwbhfjchiVVhiVVhWd2rbV0DDCy3CKbEq6dSFc7keg6mmNxkqNTwpVX4ArVYQN16RLVtbiDtSNBc6AS44yf11lEZG/QcggRkWkos5GvtqJwCrYvrvHz+VN9tPbY7JKLHf0WV3pd8SudSe5eH6e1x5KJSxtDhsiHUoVG3n93ptCIyRYbWVbtcExD4Y8FYwxVZVBV5uKw+sJBazyZmxIu83WIaCxKJNZEZEuSzv78lHDGcfCHqiBQh8nNcJGdUa7F+Cq0iU9EJt1kBMHvNcZcBqwFPmKt7SjUyBhzJXAlwPz58yfhaUVESlvuprljGvLPff1MP18/05/NVRztStI3PHLeMdDcneTRplTRD4BjG1w88o7Ux8KaX/TSOWCzWS8aQ4bD57hYk17ikbR2zPpkt2OYHzbMDxdfJxwbzBQXyU0J15VKCdf6Mk0vxhlO5E8ne3x+PKFakulMF6mZ5Jw1ysFqjKv4+mwRkUJ2GgQbY+4D6gucugb4PnA9YNNfvwa8vVA/1tqbgZsBVq9ebQu1ERGRyVUsV/FNZ42kh+tLl4geypmkPazexfNtCV5qT3L/xjjdQ3DJge5sEFz/1R58bvLWJJ+2yMXZy1LB6NaeZMGNfCGf4cA6FwcWWQaStJatY1LCJYnEWojGWohsg7aeIinhArWjUsKNfHXKQppNFpE8Ow2CrbVrdqUjY8wPgTsnPCIREdmnyj2G5dX5QelXz8jPodw1YBlIp1lLWsu7V3vSm+eSPNWS5I8vxgEvZy/z0Dtkqf9aDx4HGkImO5v8xoPdnLPMw1DC8u+2JPMLbORzjGFO0DAn6HD0qNntjEzQPpISzhLp6iESixFpf4XIxgQDw6M28bk9+MI1JIN1uDIzyqGakbXJwRocT/GUeCIy80w0O8Qca21L+scLgecmPiQREZluwn5DmFSw6hjDdafkB8nW2ryZ5O+c7c8WGYl2WR6Oxjm2IRVov9KR5PD/ThUaKfeMzCZ//Dgfpy9x09FveWzLyKa+oC9/BjcTtC+vLjzWTHXAsSnh2ojE2ohufp6W7jh21P9J+iuCmFAdBEbSwOVmvXBVzFJKOJEZZKJrgr9sjDmM1HKITcC7JjogERHZ/xhj8KU/USq8hvccVbzQRn3A4X9fX5ZXZCTaZYmnU0s80ZLgrFtGCo2EfamKfD8418/x891s7Ejy983xVP7kcKoin9+dXw1wvE2FAEMJy5YxKeEGiHRtJhKLsDmSoGMwfzbZcbnwhzIp4UYCZHewZmTZha98T19CEdnHJhQEW2vfMlkDERGR0hD2Gy5eVXwj21HzXPz9beV5RUaiMUvYnwp0/7Y5ztvuGMh7TG254f/eWs5BdS4ebYrz1835OZTnBg0e10ig7HUZFlUaFlUWTwnXNVgoJdwOIl07iGx5kS3r4iRG5YTzlpXjDtVicwqM5K1RDlQrJZzINKEUaSIiMq2EfIYT5hf/eLr0IA/HN7rG5FCuD6SC3Ac2Jbj6/vxCI46B5g8HmB1w+P26YR6O5OdQbgwZ6gMmuz7ZGMMsP8zyuzhkduGgNZG0tIzZxDdMNLaFSGwLkRbLjr5RKeGMwR+qhGBddhNfXpAcrMXxB7SJT2QfUBAsIiL7FZ/bsKzaxbIia4KvOsHHfx3pzQuQo7FUtgqAJ1sSfH/tEP05SSbcDgxckyo08rVHBnmqNZ1DOR0gL5jljAmGXY6hIWRoCDkc11h4LD1DNjuTHI1lll50p9LCbdtA9OU4Q/FRKeG8PjyhGpLB2Tlp4GrzKvIZt1LCiUyUgmAREZlxRlKxjT33+VP9XH+Kjx39I5vncguNtPakNvJtiVkySSYWhA2bPpgqNPKu/9efKjQSdpifDpQPqHY4ccHYj9SA17Cy1sXK2sLjTFpLW6/NTwfXZYnEthKNbSWyEbZ2x8c8zh8MY4J1mEBmFjk/NZxTHtZssshOKAgWEZGSY4yhutxQXc6YCnhfOcPPV87wZ3MWR2OWvuGR2dqgz9AzBPduiNPSnarId8J8F39/W+oj9YSf9NI1aLNFRhrDDqvnujhraer8UMLidY1k2pgdMMwOwJHzCi+7GIhbmnI28aXWSvcRiW0k0rGJyKYE3UOjU8K58YXyU8K50mnhsssuvP6CzydSKhQEi4iIFDCSszj/eG4O5eF0Rb6BnCUNJ853sW57agnG2mZLW5/l0oPc2SB49le78ThmZE1yyHDmUjevWZ5a4hDpSjInMLKRz+82LK0yLK0qvokvd1Z7ZFZ5O5HYdqLRdTTHhhm1hw9feQBXqBYbnJ2X4WIkJVylNvHJjKYgWEREZA95XIYFs/KXHXxxTf4Ma/+wza4/TlrLR471Zdcrb+hI8uCmJOUew2uWe+getCy4qQfHQH0gVYK6MWS47FAPr1nuYTBueWZrar3y7ECqIt94s9oZwwk/zd35QXI0NkikK0qku4lIU5LO/vxNfI7j4AtVYdOzyYWyXTi+isl7MUX2MQXBIiIie1GZx1CW3sfmGMOnXj22Ml0m1Zpj4ObX+LNp4aJdSZ7ZmqS5O3V+Q0eSo3+UKjTicWBeKDWb/KlX+zhjiZv2viSPRBPZDX1VZakgOROsL5hVvNhH10AmZ3NuSrguIl0dRFpepumFeDaXc4bXXzYmJVxesByoxrgUasj0pDtTRERkimU25VV4De88onihkYaQwx/fUJafHi6WJDMXvbY5yXm39mfbl7lThUZ+cl6q0MjL7Qke3DQSJDeGHULpinxhvyHsd3FQXfGUcFt7Ry+5SBCNNROJNRPZCtt78zfxGWPwB2flpYTLC5KDNThlIW3ikymhIFhERGQ/EfIZXntA8fRox8938egVFdm0cNF0arbq8kyhkQRX3plfaCTkg0feXsGBdS4ejsT5v1dGcijPD6dSwJV5DC7HMDdomBt0OKah8PP3DRdKCddDpCtGZPtGIuvjDI5KCef2ePGGa0gG6nCH6rJp4HKDZeMu/g8DkT2lIFhERGSGCHgNR81zcVSRTBNvPczD6Uvco3IoW+YEU8skHook+MyDg2Met+2jAWorHG59bpiHIvG8IiONYYcF4dSyi3KP4YAaFwfUFB6ftZbtfYVSwm0jGttGZBO0xAqkhAuEcIJ12OConMnp9HBORRhjii/1EClEQbCIiEiJcDuG+eHUhrtCPnGCjw8c46Upp8hIU8xSk55J/ndbglueHaYzZzLZ64L+a4IY4At/H+Sp1kReeriFs1Ip4iC1PKK2wlBbAUfMLRyoD8YtW7pHp4TrJxJ7hUjnZiKbE3QMjkoJ53LjC1WRDM7OFhcZvUbZ8ZZN/AWUGUVBsIiIiGSNl5LtulP8XHeKn+5Bm11u0TFgcdJremODlme3JvnTy3H6hlOPWVxp2PD+VJ65y29PFxoJjZSsXlnjcNrikXDE5zYsrjQsriyeEq5zgJzlFpmAeQeR2A4iTS/QHItnNxtm+y2rSKeEq8vmSs4rWR2oUkq4EqMgWERERHZL0GdYVeti1ahKeDeu8XPjmlSg2jEA0a4kvTmFRuoDhlc64Z9NCX4bizOchJMWuLJB8Oqbe4gNkpdD+egGVzaHcu+QpcJrqCyDyjLXmFLWGfGkpaV79Ca+IaKxJiLdW4g0J+noy08JZxwHf6gKArXZTXx5QXI6JZw28c0cCoJFRERkUhljqCqDqrL8IPXGnBzKSWvZ1mvpHx45f84yNy+2p2aY79sYp6XH8saDUzmSrbXUf60bA6n8yelA+Zxlbi5YkTq/fkcyu5GvMZxajnF8kTHmzmZn1ydnUsJtXU/Ty3GGRm3i8/j8eEIjm/hSwXG6Cl+oDnewGuMqvnFRphcFwSIiIrLPOcZQH8ifVb3ulPxCI/GkzS6rSFr4zKt92bRw0a4kT7bEqS03XLDCQ9cgLP9OKodyddlIRb53HO7h/BUeBuKWx7ek0sPNC5qis9kZmSA9fxNfkmislUislcgG2NZTICVcIJxOCTcqHVw6aFZKuOlDQbCIiIhMS27HEErXFnE5ho8dP7bQiLWp2VqPAz+/wJ9NCxeNWTZ1JtnRnzr/cnuSV/+sDwBDamlGY9jwuZP9nLXUTVtvMp1DORU81wcM9QGH+gBFs230D1uaYrkp4SyRrl4iXRuI7HiFyMYE3cP5m/jcbk82JVx+zuTa7Myy4/EXfD6ZXAqCRUREZL+VmVWt8BouO7R4PuGFsxzueXN5Xg7lSCyJPx0JrW1OcMlvRwqNuB2YFzTc8royjp/v5oXtCe7bmMhmvWgMGWrKDcuqXSyrLvyc1lra+23OkotMJb42orE2IhFDc2wYm7/qAl9FEFdw7Ca+zNpkV6BSKeEmgYJgERERmfGCPsMZS4qHPSctdPP0uyryioxEY5a6ipFCI++7O7/QiN8NT1xZwapaFw9uinP/xnhe/uTGkENNuUNNORw+p/Bs8nDCXyAl3ACRWIRILMrmSJLOwfxNfI7LhT9UhQ3U4YRyZpSDmTXKdTi+8gm+YjPfhINgY8z7gPcACeAua+3HJzwqERERkX2o3GM4tN7FofWFg9UrXuXhvAPcIwFyutDIvHShkce2JPjCQ0OMyszG9o8FqC53+OUzQ/w9ks6hnMl+ETYsq0rlUl44q/jMbteALZASrpNIVwfR5pdoWhcnPuqJvf4y3OmUcKlNezV5BUZcgSqMq7TnQif02xtjTgHOBw611g4aY+omZ1giIiIi00dmI199AI4ssEb448f7+PCxXlq6R7JONMWSVJWlZpI3dCT544txtvaOBKt+N/R9MpVD+doHB3iyJZk3i7yo0nBco5uw3xD2uziornCAnkhaWnvyl1xEYwkiXc1EYi1EWi3tvQU28YUq81LCudIBcnZtsj84ozfxTfSfAP8J3GitHQSw1m6b+JBERERE9j9uZyQ1G43556492c+1J/uzFfGiXUk6Bmw2yIwnYXNXkociSTrSqy6WVTm89L4AAG/4bR8bO5LZALkxZDiwzsVZS924HMPcIMwLuTl21PNm9A4VSgnXTaSri0jbRqLr4wyOSgnn9nrxhmrTKeFGguRMWjh3sAbjLr4Oe7qbaBC8HDjRGHMDMAB81Fr7eKGGxpgrgSsB5s+fP8GnFREREdn/FKuI9/lT/Xz+1NT3PUOpWeSeoZHzy6ocOgYs/25Lcs/6OL3DcOqiVBAMcPD3e+kesnnLLY5rdHHBilTe4sEEHFDtsKKm8MyutZa2vtEp4SzR2FYisa1EXoHW7viYx/kDYcyolHC5WS+c8vC03cS30yDYGHMfUF/g1DXpx1cBxwBHArcZYxZbO3qfI1hrbwZuBli9evWY8yIiIiICAa9hRU3+0ofrTx1Jm5YpHZ1bje8/DvTw8o5U5ovHtyT4w7o423o92UIiDV/vxgINoZGNe+ctd3PRqtT557alZpmPmGNYPbfwsovBuKUplllukVl60UekayORzk1ENiXoGMpPCedyu/GFqnnnv4/hrj/ePmmv0WTYaRBsrV1T7Jwx5j+B36eD3seMMUmgBmibvCGKiIiISIYxmdLRI7O6nz4pP4eytZaB9MRtwsKNa3x5WS8eeCXOolmGi1Z52NFvOeQHqUIjAS/Z2eR3H+HlwpUe+oYtj0RH0sMtqSocPuaWy85PCdeOPzj9lk1MdDnE7cApwAPGmOWAF9g+0UGJiIiIyJ4zxlCWruDsdgzvP3psoZGMMo/hNxeXZTNeZGZ6MzPNL7UnOf2Xfdn2VWWGxpDhi6f5OHuZh9aeJPdtjNMYcpgfdlhZ64zNsvGhr0/67zhREw2CfwL8xBjzHDAEvLXQUggRERERmZ7KPYZLDvQUPb+0yuHBt5aPyaEc9KVmotc2J3jLH0ZyKBtgdsDw29enCo08vy1Bbdt26sINe/tX2S0TCoKttUPAmydpLCIiIiIyzQS8hpMWFg8Z1yx2s+49FTk5lFNf56RzKD8USXDhvhrsbijtLMkiIiIiMiF+d2oj34qawuevPMIDNUVqS0+haRMEDw8P09TUxMDAwM4b76mGN0LdeXuv/xnL4u/aSMOTX8Iz1DnVgxEREZH9iDEGpmHRjWkTBDc1NREMBlm4cOHeq07SuRn6duydvmcway3tvVU08QkW/fPqqR6OiIiIyIRNm+zFAwMDVFdXz+jyfPsrYwzVFW4GwouneigiIiIik2LaBMGAAuBpLHVtdH1ERERkZphWQbCIiIiIyL4wbdYEj7bwqrsmtb9NN5477vlTLr6Sq957OWeefFz22E0/vIUXN2zm+zd+kp/f9v/4/Dd/BMCnPnAFb73ktbv83Lf8/k986Xs/w1oIVpTz/S9+kkMPXE50SyuXfeAzbN3ejjGGK9/0Oj5wxRsBeOb5l3j3VTfQ09fPwoY53PKdGwgFA3vwm6e0bG3jrR/8LPf++nt73IeIiIjITKGZ4LRLLziTW++4J+/YrXfcw6UXnMmOji4+942befTOX/DYXb/kc9+4mY7O2C73vahxHn/97Y949v7b+PQH38mVn/g8AG63i6999kP8+8Hf8c//93O++7Pb+PdLGwG44mPXceMn38+z99/GhWefwle+/4sJ/X5/fvARzjzp2An1ISIiIjJTKAhOu/jcNdx1/0MMDQ0DsCnaTPPW7Zx49Ku456//4PQTj6aqMkzlrBCnn3g0f37wkV3u+7gjD6VyVgiAY151ME0tWwGYM7uWVx28EoBgoIKVyxaxpXUbAC9tjPDqY14FwOknHsPv/nT/mH4ffGQtJ110Bee/7UMsPva1XPWFb3HL7//EUee+hYNPu4QNm6LZtn9+4BHOPvV4Wra28erXvYPDTn8DB536ev7+6JN78GqJiIiI7N8UBKdVVYY56rADufuBh4HULPAlrz0dYwxbWrfROLc+27ZhzuxssJrrM1/5Pn+896/jPs+Pb72ds085fszxTdFmnnruRY4+/CAADly+mDvueRCA/73zPqLNWwv298y/X+IHN36SdQ/+jl/+7i5e2riZx+76JVdcegHf/smtACQSCV7csJlVyxfzP3/4M2eedCxP/+VWnvnLrRx24AE7f3FEREREZhgFwTkuveCs7JKIzFKI3XHdx/6T8844qej5Bx5+nB//+na+9Mn35x3v6e3jond+lJs+95Hsut+ffP2zfO/n/8sRZ72R7t5evJ7CNb2PPPRA5syuxefzsmRBA2eklzwcvGIpm5paAHj0yeeywfWRh63ip7f9kWu/9gOeXbeeYKBit35HERERkZlAQXCO8888mfsfeownn11HX/8ARxyyCoB59XVEm1uz7ZpatjKvvm63+v7Xv1/iio9dzx0/+QbVVbOyx4eHh7nonR/lTReew+vOOS17fMXSRdz76+/xxJ//h0vPP4slCxsK9uvzjgTHjuPg83qz38fjcQDufuBhzjolteHv1cccwd9+92Pm1ddx+Yc+yy/+987d+j1EREREZgIFwTkCFeWcctxq3v7hz+XNAp950rHc+7d/0tEZo6Mzxr1/++dubTKLbGnhde/8KL/85vUsX7Ige9xayzs+ch0rly7iw+96c95jtm1PVbZLJpN8/ps/4t1vuWiPf6/7H3qMNSceDcDmpmZm11bxzje9jiveeAFPPrtuj/sVERER2V9N2xRpO0tptrdcesFZXPiOj3Dr97+YPVZVGebTH7yCI89NBaqf+dA7qaoMj3nsZ77yfVYfumrMkojrvvFD2ju6+K9Ppvp0u12svfsWHn78aX75u7s4eOVSDjv9DQB84ar3cs5pJ/Dr2//Md392GwCvO+dU3vYf5+/R79PW3oHf580ue3jwkSf4yg9+gcftJlBRxi++ef0e9SsiIiKyPzPW2n3+pKtXr7Zr167NO7Zu3TpWrly5d5+4czP07di7zzHN/Op3d9HUso2r3vu2Cfe1bvM2Vt5zySSMSkRERErKh56HcOGlnXubMeYJa+3q0cen7UywTI43XzQ1M+oiIiIi05nWBIuIiIhIyVEQLCIiIiIlZ0LLIYwxvwEy1RZmAZ3W2sMmOCYRERERkb1qQkGwtfY/Mt8bY74GdE14RCIiIiIie9mkbIwzxhjgEuDUyehPRERERGRvmqzsECcCW621LxdrYIy5ErgSYP78+Tvv8dqxeXgn5NrxJ6lPufhKrnrv5Zx58nHZYzf98BZe3LCZd735Iv7z6i8Q6+nF5XK45n3v4D/O3/WSyu07Orn4yo/z+DPPc/klr+U7N1yVPTc0NMx7P3UjDz7yBI7jcMMn3sNF557G4OAQl33g0zzx7DqqK2fxm+/fyMLGueP2tafOfvN7+eGXP03D3NkT7ktERERkf7DTINgYcx9QX+DUNdbaO9LfXwr8erx+rLU3AzdDKk/wbo5zr7v0gjO59Y578oLgW++4hy9/6gOUl/n5xTevZ9ni+TS3tnHE2W/izJOPY1Y4uEt9+/0+rv/4f/LcCxt47sX1eedu+NaPqKuu4qWHbieZTLKjMxWs//jXt1MZDrH+4T9y6x338IkbvslvfvClcfvaE/39A7R3dCkAFhERkZKy0+wQ1to11tqDCvy5A8AY4wZeB/xmbw92b7r43DXcdf9DDA0NA7Ap2kzz1u2cePSrWL5kAcsWp2av59bXUlddSVt7xy73XVFexglHHY7f5x1z7ie3/pGr3/d2ABzHoaaqEoA77n2Qt77+Nemxncb9Dz2OtXbcvnItPPpcrv7itzns9Dew+uw38eSz6zjzjf/FkuPO4we/+G223YP/eIKTjz0CgKu+8C1WnXwRh6y5hI9e941d/v1ERERE9jeTsRxiDfCCtbZpEvqaMlWVYY467EDufuBhzj/zZG694x4uee3ppJY7j3jsqecYGh5mycKxVU8yweW7L7t4l56zs6sbgE9/+Xs8+I8nWLKgge/c8Alm11azpbWNxrmpCXi32004FKC9ozMbJO+K+XPrefovt/Khz36Vyz/0WR6+/acMDA5x0Kmvz47x7gce5oIzT6Z9Ryd/uPsBXvjb7zHGZMcmIiIiMhNNRp7gN7CTpRD7i0svOItb77gHSC2FuPSC/HW/LVvbeMv7P81Pv34tjjP2pXv3ZRfvcgAMEE/EaWrZynGrD+XJe/6HY484ZFJnYM874yQADl65lKMPP5hgoILa6kp8Xm82yH348ac54ajDCYcC+H1e3vGRz/H7P91PeZl/0sYhIiIiMt1MOAi21l5urf3BZAxmqp1/5snc/9BjPPnsOvr6BzjikFXZc7HuHs697APc8In3cMwRh0zK81VXzqK8zM/rzkkl1Xj9a9bw5HMvADCvvpZocysA8XicrlgP1ZWzdqt/n88DgGMcfF5P9rjjGOKJOBs3N9E4tx6v14Pb7eaxu37Jxeeu4c77/s5Zb3rPJPyGIiIiItOTKsblCFSUc8pxq3n7hz+XNws8NDTMhe/4CJddfC4Xv2bNpD2fMYbXnv5qHnxkLQD3P/QYq5YtBlKzuD//3zsB+O1d93Pq8UeOWZoxUXc/8DBnpTcC9vT20dXdwzmnncA3rv0Iz/y7aKIPERERkf3eZKVIm3w7SWm2R4wLnPF/5UsvPJcL3/5Bbv3BV7Jtb7vzz/zt0ado74zxs3Rg+rObruewg1bkPfYHP78NgHe/9ZIx/S488ixiPT0MDQ1z+z0Pcu+v/5tVByzhS5/6MG953yf54LVfo7a6kp9+43pw3LzjjRfzlvd9kqXHn0/VrDC3/uDL2fEU62vULwvGnXqM4wLj5PzuqXN/fvAffPuGq8Fx0903yPmXf5CBwUGstXz9cx8d+1oZF5RX78ILLSIiIpLDTL95V2Ptvs9Wtnr1art27dq8Y+vWrWPlypX7fCylanBwkOOPP57R12E8ukYiIiKyvzHGPGGtXT36+PQLy2Wf8Pl8uxUAi4iIiMwkCoJFREREpORMqyB4KpZmyK7RtREREZGZZNoEwX6/n/b2dgVb05C1lvb2dvx+5Q4WERGRmWHaZIdoaGigqamJtra2qR6KFOD3+2loGFslT0RERGR/NG2CYI/Hw6JFi6Z6GCIiIiJSAqbNcggRERERkX1FQbCIiIiIlBwFwSIiIiJScqakYpwxpg3YvM+fePfVANunehCy1+j6zly6tjOXru3Mpus7c03ltV1gra0dfXBKguD9hTFmbaEyezIz6PrOXLq2M5eu7cym6ztzTcdrq+UQIiIiIlJyFASLiIiISMlREDy+m6d6ALJX6frOXLq2M5eu7cym6ztzTbtrqzXBIiIiIlJyNBMsIiIiIiVHQbCIiIiIlBwFwUUYY84yxrxojFlvjLlqqscje84Y8xNjzDZjzHM5x6qMMX8xxryc/lo5lWOUPWOMaTTGPGCM+bcx5nljzAfSx3V9ZwBjjN8Y85gx5pn09f1c+vgiY8yj6ffn3xhjvFM9VtkzxhiXMeYpY8yd6Z91bWcIY8wmY8yzxpinjTFr08em1XuzguACjDEu4LvA2cAq4FJjzKqpHZVMwM+As0Yduwq431q7DLg//bPsf+LAR6y1q4BjgPek/67q+s4Mg8Cp1tpDgcOAs4wxxwBfAr5hrV0KdADvmLohygR9AFiX87Ou7cxyirX2sJz8wNPqvVlBcGFHAeuttRuttUPArcD5Uzwm2UPW2r8BO0YdPh/4efr7nwMX7MsxyeSw1rZYa59Mf99N6sN0Hrq+M4JN6Un/6En/scCpwG/Tx3V991PGmAbgXOBH6Z8NurYz3bR6b1YQXNg8IJrzc1P6mMwcs621LenvW4HZUzkYmThjzELgcOBRdH1njPR/lz8NbAP+AmwAOq218XQTvT/vv24CPg4k0z9Xo2s7k1jgXmPME8aYK9PHptV7s3sqn1xkOrDWWmOMcgXux4wxAeB3wAettbHUhFKKru/+zVqbAA4zxswC/gCsmNoRyWQwxrwG2GatfcIYc/IUD0f2jhOstVuMMXXAX4wxL+SenA7vzZoJLmwL0Jjzc0P6mMwcW40xcwDSX7dN8XhkDxljPKQC4Fustb9PH9b1nWGstZ3AA8CxwCxjTGYSR+/P+6fjgfOMMZtILTk8FfgmurYzhrV2S/rrNlL/gD2KafberCC4sMeBZeldql7gDcAfp3hMMrn+CLw1/f1bgTumcCyyh9JrCH8MrLPWfj3nlK7vDGCMqU3PAGOMKQNOJ7Xu+wHg4nQzXd/9kLX2amttg7V2IanP2P+z1r4JXdsZwRhTYYwJZr4HzgCeY5q9N6tiXBHGmHNIrVdyAT+x1t4wtSOSPWWM+TVwMlADbAU+C9wO3AbMBzYDl1hrR2+ek2nOGHMC8HfgWUbWFX6S1LpgXd/9nDHmEFKbZ1ykJm1us9ZeZ4xZTGr2sAp4CniztXZw6kYqE5FeDvFRa+1rdG1nhvR1/EP6RzfwP9baG4wx1Uyj92YFwSIiIiJScrQcQkRERERKjoJgERERESk5CoJFREREpOQoCBYRERGRkqMgWERERERKjoJgERERESk5CoJFREREpOT8f9ycgyg9heKCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "optimized_model = layered_earth_model(v=res.x[0:2], z=[0., res.x[2]], theta=[0., res.x[3]], td= [res.x[4], res.x[5], res.x[6]])\n",
    "plt.figure(figsize=(12,4))\n",
    "optimized_model.plot_model(model_span)\n",
    "optimized_model.plot_layout(sources, receivers)\n",
    "mod.plot_model(model_span, filled=False, linestyle='--')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Edit Metadata",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}